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Introduction 


In  many  volcanic  studies,  estimates  must  be  made  of  the  changes  that  magma  and 
its  associated  gases  experience  when  traveling  through  an  eruptive  conduit  to  the  surface. 
Exsolution  of  magmatic  gas,  acceleration,  changes  in  pressure  and  temperature,  depth  of 
fragmentation,  and  final  exit  velocities  affect  such  features  as  lava  fountain  heights,  the 
ability  of  a  volcanic  column  to  convect  or  collapse,  and  the  degree  to  which  water  can 
enter  the  conduit  during  eruptive  activity.  Most  of  these  quantities  cannot  be  easily 
estimated  without  some  sort  of  numerical  model. 

This  report  presents  a  model  that  calculates  flow  properties  (pressure,  vesicularity, 
and  some  35  other  parameters)  as  a  function  of  vertical  position  within  a  volcanic  conduit 
during  a  steady-state  eruption.  The  model  idealizes  the  magma-gas  mixture  as  a  single 
homogeneous  fluid  and  calculates  gas  exsolution  under  the  assumption  of  equilibrium 
conditions.  These  are  the  same  assumptions  on  which  classic  conduit  models  (e.g., 
Wilson  and  Head,  1981)  have  been  based.  They  are  most  appropriate  when  applied  to 
eruptions  of  rapidly  ascending  magma  (basaltic  lava-fountain  eruptions,  and  Plinian  or 
sub-Plinian  eruptions  of  intermediate  or  silicic  magmas)  that  contains  abundant 
nucleation  sites  (microlites,  for  example)  for  bubble  growth. 

The  numerical  parts  of  the  program  were  written  in  Fortran  90  and  can  be  compiled 
on  any  platform  (DOS,  Unix,  Macintosh  etc.)  that  has  a  Fortran  90  compiler.  The  source 
code  to  this  model  (with  the  exception  of  certain  subroutines  taken  from  Press  et  al., 

1992)  is  posted  on  the  USGS  Cascades  Volcano  Observatory  web  site 
(littp://vulcan.wr.usgs.gov).  The  executable  version  that  is  distributed  for  the  Microsoft 
Windows®  operating  system  includes  a  graphical  user  interface  with  utilities  that 
calculate  physical  properties  of  melts,  gases,  and  melt-gas  mixtures.  Scientists  or 
educators  who  are  not  directly  interested  in  conduit  modeling  may  still  find  these  utilities 
useful.  The  program  is  free  of  charge. 


M  odel  Overview 

In  any  vigorous  magmatic  eruption,  magma  is  driven  up  a  conduit  from  some  deep 
location  to  the  Earth's  surface.  As  it  rises,  gases  come  out  of  solution,  forming  bubbles 
that  expand  to  the  point  where  they  break  the  magma  into  tiny  fragments.  Those 
fragments  become  entrained  in  a  jet  of  accelerating  gas  that  vents  violently  into  the 
atmosphere.  At  any  given  depth  in  the  conduit,  the  pressure,  velocity,  volume  fraction  of 
entrained  gas,  temperature,  and  other  characteristics  depend  on  two  sets  of  factors:  (1) 
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the  initial  pressure,  temperature,  and  composition  of  the  mixture;  and  (2)  the  length  and 
geometry  of  the  conduit. 

The  model  presented  in  this  paper  calculates  flow  properties  using  one  of  two 
methods.  Under  option  1  (Fig.  1,  left  side),  the  user  specifies  the  conduit  diameter  at  the 
base  and  top;  the  program  then  solves  for  the  pressure  and  other  flow  properties  as  a 
function  of  depth.  Under  option  2  (Fig.  1,  right  side),  the  user  specifies  the  initial 
conduit  diameter  and  a  pressure  gradient  in  the  conduit;  the  program  then  calculates  the 
conduit  geometry  required  to  produce  that  pressure  gradient. 

Option  1  Option  2 

specified  conduit  diameter  specified  pressure  gradient 

program  calculates  pressure  profile  program  calculates  conduit  geometry 


Figure  1 :  Illustration  of  the  input  variables  required  by  the  program  Conflow,  and  the  two  options 
available  for  calculating  flow  properties  as  a  function  of  depth. 

In  option  1,  the  erupting  mixture  must  satisfy  one  of  two  conditions:  (1)  if  the  exit 
velocity  is  less  than  its  sonic  velocity,  the  exit  pressure  must  equal  a  specified  final 
pressure  (usually  1  atm  at  the  exit).  Alternatively,  (2)  the  exit  velocity  must  equal  the 
sonic  velocity.  The  latter  boundary  condition  results  from  the  fact  that,  in  a  conduit  of 
constant  cross-sectional  area,  the  velocity  of  the  mixture  can  never  exceed  its  sonic 
velocity.  This  is  a  basic  tenet  of  compressible  fluid  dynamics  and  is  explained  in  a 
number  of  texts  (e.g.,  Saad,  1985).  Thus  if  the  input  pressure  at  the  base  of  the  conduit  is 
raised  above  a  certain  threshold  value,  the  erupting  mixture  will  not  be  able  to  equilibrate 
to  1  atm  pressure  by  the  time  it  reaches  the  surface.  The  exit  conditions  will  vary 
according  to  the  input  pressure,  as  shown  below: 
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Input  pressure  Exit  velocity  Exit  pressure 

<  pressure  of  static  magma  column  (pgz)  0  (no  eruption)  1  atm 

slightly  greater  than  pgz  subsonic  1  atm 

much  greater  than  pgz  sonic  >  1  atm 

The  sonic  velocity  of  mixtures  of  ash  and  gas  generally  range  from  a  few  tens  to  a  few 
hundreds  of  meters  per  second. 

In  order  to  match  the  exit  conditions  with  the  required  boundary  conditions,  the 
program  makes  successive  runs,  adjusting  the  input  velocity  after  each  one,  until  one  of 
the  two  exit  boundary  conditions  is  satisfied.  In  option  2,  successive  runs  are  not 
necessary—  an  output  pressure  of  1  atm  can  be  achieved  during  a  single  iteration  by 
calculating  the  geometry  that  gives  the  specified  pressure  gradient.  The  sonic  boundary 
condition  does  not  apply  because  the  variable  conduit  geometry  allows  the  erupting 
mixture  to  accelerate  to  supersonic  velocities. 


M  odel  Assumptions  and  Limitations 

In  constructing  the  model,  we  make  several  simplifying  assumptions.  Foremost 
among  these  is  that  the  flow  of  magma  and  exsolved  gases  is  homogeneous.  That  is,  there 
is  no  relative  movement  between  the  gas  and  liquid  as  they  ascend  the  conduit.  This 
assumption  allows  the  mixture  to  be  treated  as  a  single  fluid  whose  density,  viscosity,  and 
other  properties  are  bulk  values  for  the  mixture.  The  homogeneous-flow  assumption  is 
used  by  other  modelers  of  volcanic  eruptions,  both  mafic  and  silicic  (e.g.,  Wilson  et  al., 
1980;  Wilson  and  Head,  1981;  Head  and  Wilson,  1987;  Buresti  and  Casarosa,  1989; 
Giberti  and  Wilson,  1990),  although  its  validity  has  been  challenged  for  certain  types  of 
eruptions  or  eruptive  flow  regimes  (Vergniolle  and  Jaupart,  1986;  Dobran,  1992). 

Whether  the  gas  separates  from  the  magma  and  rises  at  a  different  velocity  depends 
largely  on  the  size  of  individual  bubbles  or  pyroclasts,  and  on  the  opportunity  for  bubbles 
to  coalesce  or  aggregate  into  larger  ones  that  rise  or  fall  more  rapidly  through  the  fluid  in 
which  they're  suspended.  The  velocity  ( u )  at  which  bubbles  rise  through  a  melt,  and 
pyroclasts  fall  through  a  gas,  can  be  calculated  from  the  following  formula  (Bird  et  al., 
1960,  p.  182): 


Mp,„  -ps) 

i  ipscD 


(1) 


where  r  is  the  bubble  or  particle  radius;  p,„  and  pg  are  melt  and  gas  densities,  respectively, 
and  CD  is  the  drag  coefficient  of  the  bubble  or  particle,  which  is  a  function  of  its  shape 
and  of  the  Reynolds  number  (Re).  For  purposes  of  this  calculation,  Re=2pur/Jl ,  where  p 
is  magma  density  (for  bubbles)  or  gas  density  (for  particles);  u  is  the  velocity  of  the 
bubble  or  particle  relative  to  the  surrounding  fluid;  and  77  is  the  viscosity  of  the 
surrounding  melt  (for  bubbles)  or  gas  (for  particles)). 

For  spheres  at  Re<~  1,  the  drag  coefficient  can  be  shown  analytically  to  be  24/Re 
(Bird  et  al.,  1960,  p.  192).  For  l</?c<~1000,  experimental  studies  have  shown  that 
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CD~IS.5/Re5li;  for  1000</?<?<~200,000,  Cd~ 0.44  for  spheres  (Bird  et  al.,  1960,  p.  192), 
and  -0.44-1.2  for  non-spherical  objects  (Hoerner,  1965;  Walker  et  al.,  1971).  The  relative 
movement  of  gas  and  melt  has  different  degrees  of  importance  below  and  above  the  depth 
at  which  magma  fragments  into  particles  that  are  entrained  in  gas.  These  differences  are 
as  follows: 

Below  the  fragmentation  depth.  Assuming  the  bubbles  to  be  spherical  and  using 
typical  values  of  the  variables  in  Eq.  (1)  for  silicic  melts,  the  ascent  velocity  of  the 

g 

bubbles  within  the  melt  is  so  small  ( <  1 0  m/s)  relative  to  ascent  velocities  of  the  melt-gas 

2  2 

mixtures  (10'M0~  m/s)  that  the  homogeneous  flow  assumption  is  reasonable.  In  basaltic 
lava-fountain  eruptions,  the  small  bubble-diameters  (0.1-1  mm  Mangan  et  al.,  1993; 
Mangan  and  Cashman,  1996)  also  produce  bubble-ascent  rates  (-10  -10’  m/s)  within  the 
melt  that  are  much  slower  than  the  ascent  rate  of  the  overall  melt  (10’  -10"  m/s).  Thus  the 
assumption  of  homogeneous  flow  should  apply  to  these  eruptions  as  well.  The  model 
tends  to  break  down  for  basaltic  eruptions  where  bubble  diameters  exceed  about  1  cm  and 
ascent  rates  are  less  than  about  10  "  m/s  (Parfitt  and  Wilson,  1995;  Vergniolle  and  Jaupart, 
1986).  Under  these  conditions  the  behavior  of  basaltic  eruptions  usually  changes  from 
fountaining  to  Strombolian  or  effusive  activity. 


-4  -3  -2  -1 

log  diameter,  meters 


Figure  2:  Terminal-fall  velocity  of  melt  particles  in  H20  gas  at  T= 900a  C,  as  a  function  of  the 
particle  diameter.  The  three  gray  regions  represent  terminal-fall  velocities  through  gases  of  three 
different  densities:  (1)  0.185  kg/m3-representative  of  H20  gas  at  900  C  and  1  atm  pressure;  (2) 

1 .849  kg/m3,  representing  the  same  gas  at  p=  1  MPa;  and  (3)  18.673,  representing  H20  gas  at 
p=10  MPa.  The  upper  boundary  of  each  gray  region  represents  the  terminal-fall  velocity  of 
spherical  clasts  of  density  2500  kg/m3.  The  lower  boundary  represents  the  terminal-fall  velocity  of 
clasts  of  density  1 000  kg/m3,  having  a  CD  of  1 .0  at  Re>  1 30. 
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Above  the  fragmentation  depth.  Figure  2  illustrates  the  terminal  velocity  of  clasts 
falling  through  gas  having  the  density  of  H2O  gas  at  10  MPa,  1  MPa,  and  0.1  MPa  (and 
900°  C),  as  a  function  of  sphere  diameter.  For  melt  particles  0.1mm  to  1  mm  in  diameter, 
their  fall  velocity  is  on  the  order  of  a  meter  per  second  or  less,  which  is  small  if  the  gas- 
ascent  velocity  is  tens  of  meters  per  second  or  more.  For  melt  particles  larger  than 
several  millimeters,  fall  velocities  could  be  meters  per  second  or  greater;  which  could  be 
a  significant  fraction  of  the  ascent  velocity  of  the  mixture. 

The  effect  of  separated  two-phase  flow  was  investigated  by  Dobran  (1992),  who 
found  that  exit  velocities  were  15-25%  higher  under  separated  flow  than  homogeneous 
flow,  and  that  exit  pressures  were  one  fourth  to  one  half  that  of  homogeneous  models. 

The  different  flow  properties  propagate  down  the  conduit  to  produce  somewhat  different 
profiles  of  velocity  and  pressure  with  depth  (Fig.  3  of  Dobran,  1992).  Because  they 
consider  the  flow  of  two  phases  explicitly,  separated-flow  models  are  more  accurate  than 
the  homogeneous  flow  models — provided  that  they  contain  appropriate  assumptions 
regarding  particle-size  distribution  and  particle  shape.  The  differences  between 
homogeneous  and  inhomogeneous  models  are  generally  smaller  than  differences  due  to 
uncertainties  in  viscosity  and  certain  other  variables  (described  later).  One  goal  of  future 
work  will  be  to  incorporate  separated  flow  into  this  model. 

Other  assumptions  are: 

1.  Gas  exsolution  maintains  equilibrium  with  pressure  in  the  conduit  up  at  least  to 
the  point  of  fragmentation.  Beyond  that  point,  the  user  has  the  option  of  shutting  off 
additional  gas  exsolution  (under  the  assumption  that  gas  exsolution  cannot  keep  pace 
with  rates  of  decompression).  This  assumption  has  been  made  in  other  models  of  conduit 
flow  (Wilson  et  al.,  1980;  Wilson  and  Head,  1981;  Giberti  and  Wilson,  1990;  Dobran, 
1992;  Papale  and  Dobran,  1993;  Papale  et  al.,  1998;  Mastin,  1995b,  1997),  though  kinetic 
calculations  and  some  experimental  work  (Mangan  and  Sisson,  1999)  suggest  that  bubble 
growth  may  be  limited  by  kinetics.  To  date,  only  the  model  by  Proussevitch  and 
Sahagian  (1996)  attempts  to  consider  the  kinetics  of  gas  exsolution  explicitly. 

2.  At  any  given  depth,  flow  properties  can  be  averaged  across  the  entire  cross- 
sectional  area  of  the  conduit.  This  assumption  simplifies  the  problem  to  a  one¬ 
dimensional  one. 

3.  The  conduit  is  vertical. 

4.  Flow  is  steady  state.  This  assumption  is  appropriate  for  eruptions  that  are 
sustained  for  many  minutes  to  hours — i.e.  basaltic  lava  fountains  and  Plinian  or  sub- 
Plinian  eruptions. 

5.  No  heat  is  transferred  across  the  conduit  walls  during  the  eruption.  For  sustained 
eruptions  through  conduits  on  the  order  of  10  m  diameter  and  1  km  long,  the  heat  flux 
through  the  conduit  walls  is  generally  two  to  five  orders  of  magnitude  less  than  that 
driven  convectively  up  the  conduit,  suggesting  that  this  assumption  is  appropriate 
(Woods,  1995). 

6.  The  gas,  melt,  and  crystals  maintain  thermal  equilibrium  during  flow.  Because 
thermal  equilibration  times  for  particles  of  the  size  typically  produced  during  volcanic 
eruptions  are  on  the  order  of  fractions  of  a  second  to  a  few  seconds  (Wilson  and  Head, 
1980),  this  assumption  is  reasonable. 
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7.  The  gas  phase  consists  only  of  H2O  gas.  With  the  exception  of  certain  alkalic 
ultramafic  magmas,  water  is  the  dominant  volatile  species  of  erupting  melts.  The 
solubility  of  water  in  melts  is  also  much  better  understood  than  that  of  other  gas  species, 
or  of  multicomponent  volatile  compositions. 

8.  There  is  no  migration  of  gas  through  the  conduit  walls.  This  assumption  limits 
applicability  of  the  model  to  cases  where  gas  generation  is  sufficiently  rapid  that  bubbles 
cannot  migrate  to  the  margin  of  the  conduit  before  they  are  released  at  the  surface.  The 
assumption  is  appropriate  for  lava-fountain  eruptions,  where  vesicle  residence  times  are 
less  than  a  minute,  and  for  silicic  high-flux  rate  eruptions,  where  the  combination  of  melt 
viscosity  and  rapid  magma  ascent  limit  the  opportunity  for  gas  to  separate  from  the  flow. 
In  slowly  fed  eruptions,  gas  escape  may  reduce  the  vesicularity  of  the  erupted  magma, 
resulting  in  the  effusion  of  lava  flows  rather  than  pyroclastic  debris  (Eichelberger  et  al., 
1986;  Jaupart  and  Allegre,  1991;  Woods  and  Koyaguchi,  1995). 

M  odel  Setup 

The  following  section  presents  the  constitutive  and  governing  equations  on  which 
the  computations  are  based. 

Governing  Equations 

Using  the  assumptions  described  earlier,  we  can  write  equations  for  conservation  of 

mass, 

(2) 

dz 


momentum. 


and  energy 


=  -pg  -pur 


/ 

R 


dp 

dz 


d  h  +  udu  +  gdz  =  0 


(3) 

(4) 


of  the  erupting  mixture.  The  variables  p,  u,  and  p  are  the  density,  velocity,  and  pressure 
of  the  mixture  in  the  conduit,  respectively;  A  is  the  conduit's  cross-sectional  area;  g  is 
gravitational  acceleration; /is  a  friction  factor  whose  value  controls  frictional  pressure 
loss  in  the  vcnrr  Bird  et  al.,  1960);  R  is  the  radius  of  the  conduit;  z  is  vertical  position 
(upwards  being  positive);  and  h  is  specific  enthalpy  of  the  magma-gas  mixture. 

Equation  2  states  simply  that  an  expansion  of  the  erupting  mixture  must  be 
accompanied  by  acceleration,  or  by  an  increase  in  cross-sectional  area  within  the  vent  in 
order  to  avoid  movement  of  material  into  a  space  already  occupied.  The  equation  is 
derived  from  the  postulate  that  the  mass  flux,  M  =puA,  is  constant  at  all  points  in  the 


“The  friction  factor  defined  by  Bird  et  al.  (1960),  used  here,  differs  by  a  factor  of  four  from  that  defined  by 
Schlichting  (1968,  p.  86)  and  used  by  Wilson  et  al.  (1980).  Therefore  the  second  term  on  the  right-hand 
side  of  Eq.  (2)  also  differs  from  the  corresponding  term  in  Eq.  (1)  of  Wilson  et  al.  (1980). 


6 


A  Numerical  Program  for  Flow  in  Eruptive  Conduits 


conduit.  Equation  3  indicates  that  acceleration  within  the  vent  may  result  from  (1) 
gravitational  forces  (first  term  on  the  right-hand  side),  (2)  frictional  forces  associated  with 
flow  (middle  term),  and  (3)  the  pressure  gradient  (right  term).  Equation  4  states  that 
changes  in  enthalpy  of  the  magma-gas  mixture  (the  first  term)  are  balanced  by  changes  in 
kinetic  energy  (the  second  term)  and  elevation  potential  energy  (the  third  term). 

By  rearranging  Eq.  (2)  as  d u=-u  ( clp/p+dA/A ),  substituting  it  into  the  term  on  the 
left  side  of  Eq.  (3),  and  rearranging,  the  following  equation  is  obtained: 


U1J  -> 

-Jr  =  Pg  +  Pu- 

dz 


/ 

R 


pu 2  dA  2  dp 

IA 

A  dz  dz 


(5) 


This  equation  can  be  made  more  tractable  by  assuming  that  the  right-hand  term,  dp/dz,  is 
approximately  equal  to  the  product  (dp/dp)s(dp/dz).  The  term  (dp/dp)s  is  the  partial 
derivative  of  density  with  respect  to  pressure  under  constant  entropy  for  the  gas-magma 
mixture.  For  homogeneous  mixtures  of  gas  dispersed  in  liquid  (or  vice  versa),  it  can 
easily  be  calculated.  Just  as  importantly,  this  quantity  is  the  squared  reciprocal  of  sound 
speed  of  the  mixture,  C  (Liepmann  and  Roshko,  1957,  p.  50).  Equation  (5)  can  therefore 
be  rewritten  as 


dp 

dz 


f 

1 

V 


.2  \ 


=  pg+pu2 


f  u2  dA 

—  -p - 

R  A  dz 


(6) 


or, 


dp 

dz 


pg+pn 2 


/ 

R 


pu2  dA 
A  dz 


1  -M2 


(7) 


where  M  is  the  Mach  number  of  the  mixture,  i.e.  its  velocity  divided  by  its  sonic  velocity. 

Equation  (7)  is  used  to  calculate  the  pressure  and  pressure  gradient  in  the  conduit. 

It  reveals  some  fundamental  properties  of  the  pressure  at  various  states  of  flow.  Under 
static  conditions,  u= 0  and  M= 0,  and  the  pressure  gradient  is  simply  -dp/dz.=pQ.  or  the 
gradient  due  to  the  static  weight  of  the  magma  column.  If  magma  is  flowing,  but  at  a 
velocity  that  is  small  relative  to  its  sonic  velocity,  M=~ 0  and  the  pressure  gradient  is  a 
function  of  the  weight  of  the  magma  column,  frictional  pressure  losses  (i.e.  the  first  and 
second  terms  in  the  numerator  on  the  right  side  of  Eq.  (7)),  and  changes  in  conduit 
geometry  (the  third  term).  As  M  approaches  1,  the  numerator  on  the  right  hand  side  of 
Eq.  (7)  must  approach  zero  in  order  to  avoid  a  singular  solution.  Setting  A=nR ,  the 
numerator  on  the  right  side  of  Eq.  (7)  must  satisfy  the  following  equality  in  order  to  be 
equal  to  zero: 


Pg  + 


fpu1 

R 


pu2  InRdR 
kR  dz 


(8) 
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Rearranging  leads  to: 


dR 

_  1 

dz 

~  2 

\ 

+  / 


(9) 


Because  the  two  terms  on  the  right  hand  side  of  Eq.  (9)  are  always  positive,  the  vent 
must  be  slightly  widening  in  the  upward  direction  in  order  for  the  sonic  velocity  to  be 
reached  (Wilson  and  Head,  1981).  In  a  constant-area  duct,  the  velocity  can  never  reach 
M=  1  regardless  of  the  driving  pressure  at  the  base  of  the  conduit  (though  from 
computational  experience  it  can  come  extremely  close).  An  increase  in  pressure  at  the 
base  of  the  conduit  will  result  in  an  increase  in  pressure  at  the  conduit  exit  and  an 
increase  in  mass  flux  (due  to  greater  density  of  the  mixture  at  the  exit).  It  will  not, 
however,  result  in  an  increase  in  the  exit  Mach  number  beyond  M=  1.  The  escaping 
magma-gas  mixture  will  equilibrate  with  atmospheric  pressure  abruptly  above  the  exit, 
through  a  series  of  shock  waves  (Liepmann  and  Roshko,  1957;  Kieffer,  1989). 

In  a  gradually  flaring  conduit,  If  M<  1  at  the  point  where  dR/dz  satisfies  Eq.  (9),  and 
the  conduit  continues  to  diverge,  the  mixture  will  decelerate  with  increasing  z  and  the 
pressure  drop  will  be  relatively  modest.  If,  on  the  other  hand,  M=  1  is  achieved  in  this 
critical  section  and  the  conduit  continues  to  diverge,  then  the  fluid  will  accelerate  to 
supersonic  velocity  and  the  pressure  will  drop  significantly  with  increasing  z.  At  this 
stage,  depending  on  the  conduit  geometry,  the  pressure  can  drop  below  p=  1  atm  prior  to 
reaching  the  conduit  exit.  If  this  is  the  case,  a  stationary  shock  wave  will  develop  within 
the  diverging  section  of  the  conduit,  through  which  the  velocity  of  the  erupting  mixture 
will  drop  abruptly  to  a  subsonic  value  and  pressure  will  rise  to  a  value  that  allows  the 
mixture  to  reach  1  atm  at  the  conduit  exit  (Saad,  1985,  p.  158). 

In  a  vent  containing  a  constant  pressure  gradient,  Eq.  (7)  is  rearranged  to  isolate  the 
variable  clAJdz  as  follows: 


dA 

A 

dp 

dz 

pu 2 

dz 

(1  -M2)  +  pg  + 


fpu 2 

R 


(10) 


This  equation  is  used  to  calculate  changes  in  cross-sectional  area  for  model  runs  in  which 
the  pressure  gradient  is  specified. 

Constitutive  Relationships 

The  following  constitutive  relationships  are  used  to  evaluate  the  terms  on  the  right- 
hand  side  of  equations  7  and  10. 

M  elt  properties 

Gas  solubdity.  The  mass  fraction  of  dissolved  gas  ( mw )  in  the  melt  is  calculated  as 
follows:  (1)  the  chemical  potential  of  water  in  the  melt  ( juw)  is  calculated  using  methods 

of  Ghiorso  and  Sack  (1995;  the  “MELTS”  method)  for  a  given  pressure,  temperature,  and 
melt  chemistry  (including  assumed  mass  fraction  dissolved  water).  (2)  The  chemical 
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potential  of  the  FLO  gas  phase  ( /!„ )  is  calculated  using  thermodynamic  relations  of  Haar 

et  al.  (1984).  (3)  The  mass  fraction  dissolved  water  ( mw )  in  the  melt  is  adjusted  until  its 

chemical  potential  equals  that  of  the  gas  phase.  The  method  of  Ghiorso  and  Sack  (1995) 
is  summarized  in  Appendix  A.  Gas  solubilities  predicted  using  this  method  are 
reasonable  approximations  to  experimental  data  (Fig.  3)  and  can  be  made  without  a  priori 
knowledge  of  solubility  for  a  given  magma  type. 


0  100  200  300 


Pressure  (MPa) 


Figure  3:  Water  solubility  (wt%) 
versus  pressure  (MPa)  for 
albite,  rhyolite,  and  mid-ocean 
ridge  basalt  (MORB), 
calculated  by  MELTS  (lines), 
and  measured  in  selected 
experiments  (symbols). 
Experimental  data  for  albite 
taken  from  Hamilton  and 
Oxtoby  (1 986);  for  basalt  from 
Hamilton  et  al.  (1 964)  and 
Dixon  et  al.  (1995);  and  for 
rhyolite  from  Holtz  et  al.,  (1995, 
“HPG8”  composition).  Melt 
compositions  used  in  the 
MELTS  calculations  to 
generate  these  lines  are  given 
in  Table  1 . 


Table  1.  Compositions  of  melts  used  for  sample  runs  in  this  document.  Oxide  compositions  are 
weight  percent  of  an  anhydrous  melt.  Basalt  is  a  MORB  composition  taken  from  Table  2  of  Dixon 
et  al.,  (1 995).  Pinatubo  melt  composition  taken  from  Luhr  &  Melson  (1 997).  Rhyolite  composition 
represents  a  haplogranitic  melt  (HPG8)  characterized  for  its  solubility  (Holtz  et  al.,  1995)  and 
rheologic  properties  (e.g.,  Dingwell  et  al.,1996;  Hess  and  Dingwell,  1996). 


property 

basalt 

Mt.  St. 
Helens 

Pinatubo 

rhyolite 

temperature  (C) 

1200. 

930. 

780. 

750. 

Si02* 

50.80 

74.18 

77.53 

76.69 

AI2O3 

13.70 

14.83 

12.81 

12.91 

Fe203 

0.00 

0.00 

0.00 

0.61 

FeO 

12.40 

2.10 

0.80 

0.55 

MgO 

6.67 

0.51 

0.23 

0.04 

CaO 

11.50 

2.39 

1.30 

0.29 

Ti02 

1.84 

0.37 

0.14 

0.10 

Na20 

0.68 

5.24 

4.16 

4.20 

k2”o 

0.15 

0.37 

2.98 

4.61 

Conflow  uses  the  MELTS  procedure  to  calculate  solubility  through  a  range  of 
pressure  at  the  beginning  of  the  model  run.  It  then  uses  a  least-squares  routine  (Press  et 
al.,  1992,  p.  659)  fits  the  results  to  the  equation  below: 

mw  =  (Jpp  (11) 
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where  mw  is  the  mass  fraction  dissolved  water  in  the  melt  (distinguished  from  mw,  the 

mass  fraction  water  in  the  melt+gas+crystal  mixture);  and  <7  and  /3  are  constants  whose 
values  are  determined  by  the  best-fit  procedure.  Equation  (1 1)  is  used  to  compute  mw 
except  in  cases  where  the  total  water  in  the  system  is  not  sufficient  to  saturate  the  magma. 
2)  Mass  fraction  crystals.  The  volume  fraction  crystals  in  melt  ( V  x )  is  given  as 


input  to  the  model.  The  mass  fraction  crystals  in  the  melt  ( mx ,  distinguished  from  the 
mass  fraction  crystals  in  the  mixture,  mx )  is  calculated  from  the  formula: 

PM, 


m , 


pyx  +  pj\-vx) 


(12) 


where  the  crystal  density  is  given  as  input  to  the  program,  and  melt  density  is  calculated 
using  method  of  Ghiorso  and  Sack  (1995;  see  Appendix  A)  after  computing  the  dissolved 
water  content  of  the  melt. 

M  ixture  properties 

Mass  fractions  gas,  crystals,  and  melt.  The  mass  fraction  gas  in,  the  total  system  is 
equal  to  the  total  water  in  the  system  minus  that  dissolved  in  the  melr:  It  is  calculated 
from  the  following  equation: 


m„  = 


=  m 


w 


(i 


m 


■  m 


where  mx  and  mm  are  the  mass  fractions  of  crystals  and  melt  in  the  mixture,  respectively. 
By  substituting  mx  (l  -  mg )  =  mx  and  rearranging,  we  get: 


_  mw  —  (l  -  nix  >hu, 
1  -  (l  -  mx  )mw 


(13) 


Mass  fractions  of  the  crystals  and  melt  are  then  calculated  as  follows: 

mx  =  (l  -  mg  )mx  (14) 

mm=l  -mg-mx  (15) 

Volume  fractions.  Volume  fractions  ( V)  of  the  three  phases  are  calculated  as 
follows: 


3  In  this  model,  the  water  incorporated  into  minerals  is  ignored. 
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V  = - - -  (16) 

ms  ,  »k,  |  mx 

Pg  Pm  Px 


where  the  subscript  i  refers  to  one  of  the  three  phases;  gas  (g),  melt  ( m ),  or  crystals  (x). 
The  density  of  the  gas  is  calculated  using  relations  of  Haar  et  al.  (1984)  at  the  current 
absolute  temperature  ( T )  and  p.  Densities  of  the  other  phases  are  calculated  as  explained 
for  Eq.  (12). 

Density.  The  bulk  density  of  the  mixture,  p,  is: 


P  = 


1 

(mgVg  +  mmV  m  +mxV  J 


(17) 


Pressure  (MPa)  at  saturation 
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Figure  4:  (a)  Variation  in 
viscosity  with  dissolved 
water  content  for  a  rhyolitic, 
water-saturated  melt 
(HPG8)  with  r=750°C  using 
relations  of  Hess  &  Dingwell 
(1996),  and  Shaw  (1972). 
Composition  of  this  melt  is 
given  in  Table  1. 


Friction  factor 

The  frictional  resistance  of  single-phase  fluids  flowing  in  cylindrical  conduits  is 
well  known  from  experimental  data  (e.g.,  Bird  et  al.,  1960,  p.  186).  Frictional  resistance 
is  generally  expressed  as  a  friction  factor,/,  defined  as  the  force  resisting  flow  through  a 
unit  length  of  a  conduit,  normalized  to  the  surface  area  of  the  conduit  in  that  path  length 
and  to  the  kinetic  energy  per  unit  volume  of  the  flowing  mixture  (Bird  et  al.,  1960,  p. 
181).  Following  previous  investigators  (Wilson  et  al.,  1980;  Giberti  and  Wilson,  1990; 
Dobran,  1992),  we  calculate /from  the  following  equation: 


,  16  ,  16 p 

f  -  —  +  fo~  —z  +  fo 

Re  puD 


(18) 


where  D  is  the  conduit  diameter,  p  is  the  viscosity  of  the  mixture,  and  Re  is  the  Reynolds 
number,  defined  as  puD/r).  The  variable/)  is  an  empirically  derived  factor  related  to  the 
roughness  of  the  conduit  walls.  In  Conflow  it  is  assumed  to  be  0.0025. 
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Figure  5:  Pressure  profile 
in  a  100-m-diameter 
conduit  for  high-silica 
rhyolite  (composition  given 
in  Table  3)  with  choked 
flow  at  the  exit,  p=  1 30  MPa 
at  5  km  depth.  The 
lowermost  three  curves  use 
viscosity  relations  of  Shaw 
(1972)  and  input 
temperature  of  740  C;  all 
other  lines  use  viscosity 
relations  of  Hess  and 
Dingwell  (1996)  and  input 
temperatures  as  indicated. 
The  significance  of  the 
variable  N  is  explained  on 
p.  18-19. 

Pressure  (MPa) 

For  laminar-flow  conditions  (Re<~2000),  which  characterize  nearly  the  entire 
conduit  below  the  fragmentation  depth,  the  left-hand  term  on  the  right  side  dominates  Eq. 
(18).  The  velocity  u  is  determined  with  knowledge  of  p  and  A  using  the  continuity 
equation  (Eq.  (2)),  and  D  is  calculated  or  specified.  At  Reynolds  numbers  typical  for 
turbulent  flow  (generally,  the  conduit  section  above  the  fragmentation  depth),  the  friction 
factor/is  determined  primarily  by  the  right-hand  term ,/G,  in  Eq.  (18).  Experimental 
values  of/o  range  from  about  0.001  to  0.02;  values  of  around  0.0025  are  commonly  used 
to  model  flow  in  rough-walled  eruptive  conduits  (Wilson  et  al.,  1980;  Giberti  and  Wilson, 
1990),  and  we  use  that  value  here.  Variations  in/0  between  0.002  and  0.02  have  an 
insignificant  effect  on  conduit  pressures  for  basaltic  and  silicic  magmas. 

Viscosity  of  the  melt.  The  viscosity  (t))  of  the  mixture  varies  greatly  during  ascent 
due  to  vesiculation,  fragmentation,  heating  or  cooling,  and  the  exsolution  of  dissolved 
water.  The  method  of  Shaw  (1972)  remains  the  only  one  that  allows  viscosity  to  be 
calculated  for  any  given  melt  composition  and  temperature.  Conflow  uses  this  method 
for  all  melts  containing  less  than  70%  SiCH  by  weight  (anhydrous).  The  method  of  Shaw 
(1972)  assumes  that  the  viscosity  of  the  melt  (77, „)  is  Arrhenian,  i.e.  that  it  obeys  the 
relation: 


log(77m)  =  A  +  B/T  (19) 

where  A  and  B  are  empirical  constants  which  are  functions  of  melt  composition,  and  T  is 
absolute  temperature.  The  method  for  calculating  A  and  B  is  described  in  Shaw  (1972). 
For  melts  containing  more  than  70%  silica,  Conflow  uses  a  non- Arrhenian  viscosity 
relation  published  by  Hess  and  Dingwell  (1996)  . 


4  In  Hess  and  Dingwell’ s  original  equation,  water  content  was  expressed  in  weight  percent.  I  have 
converted  the  equation  so  that  water  content  is  expressed  in  mass  fraction  of  the  melt. 
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(20) 


log(77j  =  .2911  +  0.8331n(mj  + 


- 1304- 2368  ln(mj 
T  -(344  +  32.25  ln(mj) 


where  viscosity  is  in  Pascal  seconds  and  temperature  (7)  is  in  Kelvin.  This  relation  gives 
substantially  lower  viscosities  than  Shaw  (1972)  at  dissolved  water  contents  of  1.5-2. 5 
wt%,  and  higher  viscosities  for  water  contents<0.5  wt%  (Fig.  4).  For  silicic  melts,  these 
viscosity  relations  produce  much  different  pressure  profiles  (Fig.  5). 

Effect  of  crystals  on  viscosity.  The  presence  of  crystals  in  the  melt  generally 
increases  resistance  to  flow.  Studies  in  the  engineering  literature  (e.g.,  Einstein,  1906, 
1911;  Hess,  1920;  Eilers,  1943;  Roscoe,  1952;  Gay  et  al.,  1969;  Jeffrey  and  Acrivos, 

1976)  and  in  the  Earth  Sciences  (e.g.,  Shaw,  1965,  1969;  Kerr  and  Lister,  1991;  Pinkerton 
and  Stevenson,  1992;  LeJeune  and  Richet,  1995)  suggest  that  the  rheology  of  a  crystalline 
melt  depends  on  the  volume  fraction  crystals  and  on  their  shape.  In  general,  melts  remain 
Newtonian  as  long  as  crystals  make  up  less  than  a  few  tens  of  percent  of  the  mixture  (by 
volume;  Pinkerton  and  Stevenson,  1992;  Lejeune  and  Richet,  1995).  At  higher  volume 
fraction  the  mixture  develops  a  yield  strength  (Kerr  and  Lister,  1991;  Pinkerton  and 
Stevenson,  1992).  The  Einstein-Roscoe  equation  is  generally  used  to  calculate  the 
viscosity  ( Tjm+X )  of  a  crystal-melt  mixture  (e.g..  Marsh,  1981;  Pinkerton  and  Stevenson, 
1992;  Lejeune  and  Richet,  1995): 


m+x 


( 


v 


max 


x-2.5 

J 


(21) 


where  K,U1X  is  the  volume-fraction  crystals  at  which  maximum  packing  is  achieved.  For 
roughly  equant  crystals  packed  irregularly,  Marsh  (1981)  suggests  that  Vmax~0.6.  That 
value  is  used  in  Conflow. 

Neither  Conflow  nor  any  other  current  volcanic-conduit  model  handles  the  non- 
Newtonian  rheologies  of  melts  at  high  crystal  fractions.  If  you  enter  a  crystallinity  above 
30%  while  using  Conflow,  you  will  receive  the  following  warning: 


You  entered  a  crystallinity  of  _ %. 

At  crystallinities  above  about  30%, 

the  rheological  law  used  in  this  model  becomes  inaccurate. 

Do  you  wish  to  continue?  (y/n) : 

If  you  enter  more  than  59  volume  percent  crystals,  you  will  receive  the  following  error 
message: 


You  entered  a  crystallinity  of  _ %. 

This  program  cannot  handle  crystallinities  above  59%.  Program  stopped. 

Effect  of  bubbles  on  viscosity.  Although  numerous  investigators  have  measured 
the  rheology  of  foams  and  emulsions  (see,  for  example,  issues  of  Journal  of  Rheology  or 
Journal  of  Colloid  and  Interface  Science ),  few  studies  have  explicitly  addressed  the 
rheology  of  bubbly  melts;  and  those  few  studies  have  reached  rather  variable  conclusions 
Experiments  on  GeOo  containing  from  0.8  to  5.5  vol%  air  bubbles  (Stein  and  Spera, 
1992)  found  substantial  increases  in  viscosity  with  bubble-volume  fraction;  but 
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oscillatory-strain  experiments  on  extremely  high- viscosity  rhyolite  (Bagdassarov  and 
Dingwell,  1992)  found  decreases  in  viscosity  with  increasing  bubble  volume  fraction. 

Manga  et  al.  (1998)  appear  to  offer  an  explanation  of  these  discrepancies  by 
explaining  that  the  bulk  viscosity  is  a  function  of  the  capillary  number: 

Ca  =  (22) 

X 

where  e  is  the  shear  strain  rate  of  the  suspension,  r  is  the  average  undeformed  bubble 
radius,  and  X  is  the  surface  tension  of  the  liquid.  For  Ca«  1 ,  bubbly  melts  tend  to  be 
more  viscous  than  the  liquid  alone,  while  for  Ca>~  1  they  tend  to  be  less  viscous.  In 
general,  bubbly  mafic  melts  have  lower  capillary  numbers  than  silicic  melts;  hence  mafic 
melts  should  tend  to  increase  in  viscosity  with  vesicularity;  silicic  ones  to  decrease. 

For  low  capillary  numbers,  the  amount  by  which  viscosity  changes  with  volume 
fraction  gas  is  not  well  defined.  Taylor  (1932)  suggests  that  viscosity  of  emulsions 
containing  sparse  fluid  droplets  varies  as  r\^i)m+x(\+Vg)  (for  rjg«rjm+x).  Dobran  (1992) 
uses  the  relation  Tj=Tjm+x/(  1  -  Vg)  (  for  ?)„«?)„,).  The  increase  in  viscosity  with  V„  given  by 
Dobran’ s  equation  is  more  modest  than  that  for  hard  spherical  inclusions  (Roscoe,  1952), 
but  greater  than  that  for  liquid  droplets  (Taylor,  1932),  or  for  bubbles  at  Ca= 0.3 
calculated  numerically  for  silicate  melts  (Manga  et  al.,  1998).  The  viscosity  predicted  by 
relations  of  Dobran  (1992)  is  also  significantly  less  than  that  used  by  Jaupart  and  Allegre 
(1991)  based  on  experimental  measurements  of  Sibree  (1933)  for  bubbly  liquids;  and  is 
less  than  relations  derived  by  Stein  and  Spera  (1992)  for  bubbly  GeC>2  melts.  As  pointed 
out  by  Manga  et  al.  (1998),  the  experiments  of  Sibree  (1933)  may  not  be  applicable  to 
silicate  melts  because  an  organic  colloid  produces  adsorption  layers  on  the  bubble  walls 
of  his  liquids.  Similarly,  the  results  of  Stein  and  Spera  (1992),  which  produce  viscosities 
greater  than  that  expected  for  hard  spheres,  may  have  been  affected  by  quenching  or 
crystallization  along  bubble  walls  (Manga  et  al.,  1998).  In  the  absence  of  more  definitive 
data,  the  relation  of  Dobran  appears  to  be  a  reasonable  approximation  of  77  for  Ca«l. 

For  Ca»l,  the  bulk  viscosity  approaches  a  theoretical  limit  of  T)=T)m+x(\-Vg)  (Manga 
et  al.,  1998).  Based  on  this  relationship  and  that  of  Dobran  (1992),  one  could  postulate  a 
bulk  viscosity  given  by: 


(23) 

where  N  is  an  adjustable  constant  that  varies  from  ~1  (for  Ca«  1 )  to  —1  (for  Ca»  1). 

For  basaltic  melts,  pressure  and  velocity  profiles  are  not  especially  sensitive  to  the 
particular  viscosity- vesicularity  relationship  (Mastin,  1995b,  Fig.  7).  For  silicic  melts,  the 
nature  of  the  viscosity- vesicularity  relationship  could  dramatically  affect  flow  properties, 
as  described  later. 

Conflow  calculates  bulk  viscosity  using  Eq.  (23)  and  an  estimated  value  of  N  as 
follows: 


N  =  —  tan  *(5  ■  log(Cr/)) 
K 


(24) 
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where  the  Capillary  number  is  calculated  as  described  in  Appendix  B.  This  relation  has 
been  chosen  simply  because  it  changes  gradually  from  asymptotic  values  of -1  at  Ca«  1 
to  1  at  Ca»  1.  Additional  research  may  lead  to  an  improved  understanding  of  the 
rheology  of  bubbly  liquids. 

Viscosity  above  the  fragmentation  depth.  At  high  vesicularity,  the  bubbly 
suspension  breaks  up  into  a  gas  entraining  particles  of  melt.  For  the  viscosity  of  the 
fragmented  mixture,  Conflow  uses  the  following  relation  (Dobran,  1992): 


V=VS 


1-1.56 


0.62 


(25) 


Fragmentation 

The  point  at  which  the  melt  breaks  up  into  small  fragments  entrained  by  gas  is  of 
great  importance  in  controlling  dynamics  of  conduit  flow.  That  point  has  traditionally 
been  assumed  to  take  place  when  \/g=0.75,  the  gas  volume  fraction  at  which  spherical 
bubbles  reach  a  closest-packing  structure  (Sparks,  1978).  Several  conduit  models  (e.g., 
Wilson  and  Head,  1980;  Wilson  et  al.,  1981;  Giberti  and  Wilson,  1990;  Dobran,  1992) 
use  14=0.75  as  a  criterion  for  fragmentation.  Conflow  uses  this  criterion  as  well. 

Recent  studies  have  explored  more  physically  based  fragmentation  mechanisms, 
including  the  degree  of  overpressure  in  bubbles  (Alidibirov,  1994;  Alidibirov  and 
Dingwell,  1996)  and  shock-wave  propagation  (Barmin  and  Melnik,  1993).  Papale  (1999) 
suggests  that  fragmentation  takes  place  when  the  extensional-strain  rate  (  £__ )  within  the 
conduit  exceeds  that  which  can  be  accommodated  by  viscous  flow.  Papale’ s  criterion  is 
expressed  mathematically  as  follows: 


du  k  , 
e77  =  —  >-  =  k  — 
dz  z  77 


(26) 


where  (dw/dz)  is  the  vertical  velocity  gradient,  k  is  an  empirical  constant,  z  is  the  magma 
structural  relaxation  time,  77  is  the  viscosity  of  the  mixture,  and  G„  is  the  “elastic” 
modulus  of  the  bubbly  liquid  at  infinite  frequency.  Using  values  of  fc=0.01,  Gj=  25  GPa, 
and  77=  T)m+x/(l-Vg),  Papale  tested  this  criterion  for  rhyolitic,  dacitic,  and  basaltic  conduit 
flow.  He  found  that  fragmentation  took  place  when  the  volume  fraction  gas  range  from 
about  0.62  to  0.93,  with  higher  values  for  mafic  melts  and  lower  values  for  silicic  ones. 
Using  Papale’ s  fragmentation  criterion,  Mastin  (1999)  found  that  the  depth  of 
fragmentation  and  other  flow  properties  in  the  conduit  vary  dramatically  depending  on  the 
exact  relation  for  77  (i.e.  constant  versus  variable  A;  Fig.  7).  Fragmentation  depths  and 
flow  properties  are  also  highly  sensitive  to  values  of  k  and  G,„ 


M  odel  Setup 


15 


E 


a 

o 

Q 


p,  MPa  Volume  fraction  gas  Velocity,  m/s 


Figure  6:  Comparison  of  results  for  conduit  flow  of  Mount  St.  Helens  magma  using  a 
fragmentation  criterion  of  1/3=0. 75  (solid  line)  and  a  strain-rate  based  fragmentation  criterion  (the 
three  other  lines).  The  parameters  used  to  generate  these  profiles  are  described  in  the  text.  The 
long-dashed  line  represents  flow  profiles  generated  with  the  values  of  k,  G«,  and  77  used  by 
Papale,  including  calculation  of  r)m  using  relations  of  Shaw  (1972).  The  dot-dashed  line  represents 
profiles  generated  using  the  same  relations,  but  calculating  r/m  using  relations  of  Hess  and 
Dingwell  (1 996).  The  dotted  line  uses  the  same  relations  as  the  dot-dashed  line,  but  with  a  value 
of  77  (Eq.  23)  that  depends  on  capillary  number.  The  initial  melt  composition,  pressure  and 
temperature  used  in  these  models  are  listed  in  Table  1.  Conduit  diameter  =60  m. 

With  some  minor  modifications  of  the  source  code  (described  in  Appendix  C), 
Conflow  is  capable  of  using  Papale’ s  fragmentation  criterion.  The  Papale  fragmentation 
criterion  is  not  available  in  the  standard  executable  program  because  there  is  still  a  great 
deal  of  uncertainty  regarding  both  the  appropriate  numerical  values  (or  mathematical 
relations)  for  k,  GC,  and  77,  and  their  appropriate  definitions.  For  example,  the  elastic 
modulus  G...  used  to  calculate  the  brittle  failure  of  elongating  glass  fibers  is  the 
elongational  (or  Youngs)  modulus;  but  in  eruptive  conduits  with  rigid  conduit  walls,  the 
bulk  modulus  may  be  more  appropriate.  For  glass,  the  former  ranges  from  about  25  to  78 
GPa,  while  the  latter  is  about  twice  that  (Bansal  and  Doremus,  1986).  Papale  uses  a 
constant  value  of  Gj,  but  for  a  bubbly  mixture  the  value  of  G,  must  decrease  dramatically 
as  gas  volume  fraction  increases. 
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Similarly,  for  viscosity,  Papale  (1999)  used  the  standard  Newtonian  shear  viscosity 
for  77,  though  the  criterion  for  brittle  failure  of  glass  fibers  uses  the  elongational  viscosity, 
which  relates  extensional  strain  rate  to  tensile  stress.  Elongational  viscosity  is  generally 
about  three  times  the  shear  viscosity  (Webb  and  Dingwell,  1990).  For  conduits  with  rigid 
walls,  a  third  viscosity,  termed  the  volumetric  viscosity,  may  be  the  most  important.  The 
volumetric  viscosity  relates  volumetric  changes  of  the  mixture  to  pressure  differential 
(Thomas  et  al.,  1994;  Kaminski  and  Jaupart,  1997).  Its  value  is  not  well  established. 

Conflow  users  who  employ  Papale’ s  criterion  should  do  so  after  devoting  some 
careful  thought  to  the  parameters  involved  and  their  significance.  The  high  sensitivity  of 
flow  properties  to  such  factors  as  viscosity  may  reflect  a  real  instability  in  eruptive 
dynamics;  that  is,  minor  changes  in  crystal  content,  temperature,  or  melt  chemistry  during 
an  eruption  may  produce  significant  changes  in  eruptive  dynamics.  Additional  study  of 
the  criteria  that  control  fragmentation,  using  this  model,  would  be  a  fruitful  avenue  of 
research. 


M  ach  number 

The  Mach  number  of  the  mixture  is  its  velocity  divided  by  the  mixture's 
(approximate)  sonic  velocity  (C).  The  latter  is  defined  as 


C2  = 


^  dp  x 
dp 


Js 


(27) 


where  the  subscript  s  indicates  constant  entropy  conditions.  This  equation  can  also  be 
written  in  terms  analogous  to  seismic  velocity  equations,  as 


C2 


K_ 

P 


(28) 


where  K  is  the  bulk  modulus  of  the  mixture  under  adiabatic  (constant-entropy) 
conditions.  For  a  dispersed  mixture  of  particles  in  gas,  the  bulk  modulus  is: 


K 


v 


v 


K. 


■  + 


K 


■  + 


(29) 


where  14,  Vg,  and  Vx  represent  the  volume  fraction  of  the  three  phases,  and  Km,  Kg  and  Kx 
their  bulk  moduli.  The  bulk  modulus  of  the  crystals  is  assumed  to  be  approximately  105 
MPa.  The  bulk  modulus  of  unvesiculated  magma  is  calculated  using  the  method  of 
MEETS,  given  in  Appendix  A  (for  the  melt,  we  assume  that  the  isothermal  bulk  modulus 
is  essentially  equal  to  the  isentropic  bulk  modulus).  The  bulk  modulus  of  the  gas  phase 
can  be  calculated  from  ideal  gas  relations: 


dp  ' 

„  CPS 

r  dp  ' 

l*J 

=  Pg - 

8  c 

s  vg 

&  ;T 

(30) 
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where  the  subscripts  s  and  T  refer  to  constant  entropy  and  constant  temperature, 
respectively.  All  of  the  terms  on  the  right-hand  side  can  be  calculated  using  the  Haar  et 
al.  (1984)  equation  of  state  for  H20. 

Numerical  Procedure 

For  the  case  of  specified  cross  sectional  area  in  the  conduit,  all  terms  on  the  right- 
hand  side  of  Eq.  (7)  can  be  determined  as  long  as  the  pressure  and  velocity  at  the  base  of 
the  conduit  are  specified.  By  calculating  dp/dz  from  Eq.  (7),  a  new  pressure  can  be 
extrapolated  to  a  higher  point  in  the  conduit.  The  continuity  equation,  Eq.  (2),  as  well  as 
the  constitutive  relations  in  equations  1 1-30  and  the  appendices,  can  be  used  to  evaluate 
density,  velocity,  friction  factor,  and  Mach  number  at  this  new  depth.  Using  these  values, 
a  new  dp/dz  can  be  evaluated  using  Eq.  (7),  and  the  procedure  is  repeated  to  the  top  of  the 
conduit.  For  the  case  of  constant  pressure  gradient,  the  procedure  is  the  same  except  that 
a  new  gradient  in  cross-sectional  area  is  evaluated  at  each  depth  using  Eq.  (10),  rather 
than  a  new  pressure  gradient  using  Eq.  (7). 

The  integration  is  carried  out  using  a  Cash-Carp  method  with  automatic  quality 
control  that  adjusts  the  vertical  step  size  to  concentrate  calculations  at  points  where 
properties  are  changing  most  rapidly  (Press  et  al.,  1992). 

Temperature  changes  at  each  depth  are  calculated  using  the  following  equation, 
which  is  modified  from  Eq.  (4): 

h  =  hi+^(u; -u2)+g(z1-z)  (31) 

where  ui:  zi,  and  hi  are  the  velocity,  elevation,  and  enthalpy  at  the  base  of  the  conduit; 
and  u,  z,  and  h  are  the  same  variables  at  the  current  depth.  The  initial  specific  enthalpy 
(hi)  is  calculated  for  the  known  pressure,  temperature,  and  composition  of  the  melt  using 
the  formula: 


h  =  "(A  +  mA,  +  mA  (32) 

where  hg,  hm,  and  hx  are  the  specific  enthalpies  of  the  gas,  melt,  and  crystals,  respectively. 
The  specific  enthalpies  of  the  gas  and  melt  are  calculated  using  methods  of  Haar  et  al. 
(1984)  and  of  Ghiorso  and  Sack  (1995),  respectively:  the  specific  enthalpy  of  the  crystals 
is  calculated  using  the  following  simplified  equation: 

h  =  cxT  +  —  (33) 

Px 

where  the  specific  heat  (cx)  and  density  (px)  of  the  crystals  (assumed  constant)  are  given 
as  input  to  the  program. 

At  any  depth  above  the  base  of  the  conduit,  the  elevation  and  velocity  are  used  in 
Eq.  (31)  to  calculate  a  new  enthalpy  of  the  erupting  mixture.  For  the  known  pressure  and 
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composition  at  that  depth,  the  program  adjusts  the  temperature  of  the  mixture  until  its 
enthalpy,  calculated  using  Eq.  (32),  equals  that  predicted  by  Eq.  (31). 


Pressure,  MPa 


Figure  7:  Comparison  of  certain  melt  or  flow  properties  calculated  by  Conflow  and  calculated  by 
independent  methods,  for  a  Mount  St.  Helens  magma  at  830  C,  with  4.6  wt%  total  water.  These 
plots  indicate  the  accuracy  of  Conflow  in  calculating  these  properties.  In  the  upper  plot,  the  sonic 
velocity  calculated  by  Conflow  (solid  line)  is  compared  with  that  calculated  using  the  equation: 


\-v 


g 


where  the  bulk  moduli  of  the  gas  (Kg)  and  of  the  melt+crystals  (Km+X)  are  calculated  by  Conflow. 

In  the  middle  plot,  the  gas  density  calculated  by  Conflow  is  compared  with  that  calculated  using  a 
Fortran  program  for  steam  properties  provided  by  J.S.  Gallagher  of  the  National  Bureau  of 
Standards.  For  the  sake  of  comparison,  we  also  plot  density  of  an  ideal  gas  (squares)  having  the 
molecular  weight  of  water.  In  the  lower  plot,  the  enthalpy  of  the  melt  calculated  by  Conflow  is 
compared  with  the  values  for  the  same  p,  7~,  and  composition  calculated  by  the  web-based 
MELTS  calculator  (http://weber.  u. Washington. eduZ-ghiorso/  . 


Testing  the  M  odel 

The  tests  presented  in  this  section  illustrate  two  points:  (1)  that  the  model  correctly 
calculates  various  properties  of  the  erupting  mixture  as  set  forth  in  the  constitutive 
equations;  and  (2)  that  the  model  correctly  calculates  flow  properties  for  certain  end- 
member  situations  for  which  analytical  solutions  exist. 

In  addressing  point  (1),  we  do  not  attempt  to  show  exhaustively  that  every  property 
calculated  by  Conflow  is  accurate:  however  in  Fig.  8  we  illustrate  the  accuracy  of  a  few 
key  parameters  (sound  speed,  gas  velocity,  melt  enthalpy)  by  comparing  values  calculated 
by  Conflow  with  independent  calculations.  The  results  compare  well  (as  one  would 
expect). 
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To  address  point  (2),  we  calculate  conduit  flow  under  two  end-member  conditions: 
(1)  isothermal  flow  of  a  single-phase  melt  through  a  conduit  of  constant  cross  sectional 
area;  and  (2)  flow  of  a  perfect  gas  through  a  frictionless  conduit.  The  overall  results  of 
these  end-member  tests  depend  on  each  of  the  flow  properties:  therefore  in  addressing 
point  (2)  above,  we  are  implicitly  testing  point  (1). 
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Figure  8:  Comparison  of  flow  properties  calculated  by  Conflow  (solid  lines)  and  the  analytical 
solution  (dashed  lines)  for  isothermal,  laminar  flow  of  an  incompressible  liquid  (Kilauean  basalt  at 
1145°C)  up  a  1-km  long  vertical  conduit,  5  meters  in  diameter.  In  the  left  plot,  the  results 
calculated  by  Conflow  are  indistinguishable  from  those  of  the  analytical  equation.  Conflow 
considers  changes  in  density  (middle  plot)  and  temperature  (right-hand  plot),  which  are  not 
considered  by  the  analytical  equation;  however  those  changes  do  not  significantly  affect  the 
calculations  of  pressure. 

Steady,  Isothermal  Flow  through  a  Conduit  of  Constant  Cross-sectional  Area 

The  continuity  equation  (Eq.  (2))  for  this  case  reduces  to  p=constant.  Equation  5 
reduces  to 


dp 

dz 


(34) 


Substituting /=16/i?e,  and  considering  that  Re=2pur/i],  the  equation  can  be  rewritten  as 
follows: 


dp 

dz 


This  is  easily  integrated  to  give: 


Pf  -Pi 


f  Spu  \ 

pg+-^r  l zf-z i) 
r  ) 


(35) 


(36) 
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Figure  9:  Flow  properties 
in  a  5-m  long  conduit 
containing  pure  H20  gas, 
having  a  diameter  at  the 
base  of  5  m  and  a 
pressure  at  the  base  of 
0.2  MPa.  Solid  lines  give 
the  flow  properties 
calculated  by  Conflow: 
triangles  give  the  flow 
properties  for  a  perfect 
gas  with  y=1 .249. 
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where  the  subscripts/ 


and  1  refer  to  the  final  and  initial  values,  respectively,  of  p  and  z.  Figure  9  compares  the 
pressure  profile  (left),  melt  density  (center)  and  temperature  (right)  calculated  for  a 
Hawaiian  basalt  using  Conflow,  and  using  the  analytical  solution  with  a  volatile-free 
magma,  initially  at  1 145°C  (21.40  Pa  s  viscosity),  flowing  at  0.001  m/s  through  a  2-cm- 
diameter  conduit  (the  conduit  diameter  and  velocity  had  to  be  adjusted  to  ensure  that  flow 
was  laminar).  The  pressure  profile  given  by  Conflow  (solid  line)  matches  the  analytical 
solution  (dashed  line)  closely,  but  not  exactly.  The  discrepancy  is  due  to  adiabatic 
changes  in  temperature  of  the  magma  (-0.25°  cooling  after  1000  m  of  flow  (middle  plot), 
which  increases  its  viscosity  by  about  0.08  Pa  s  (lower  plot)  and,  combined  with 
decompression  effects,  decreases  its  density. 

Choked  Flow  of  a  Frictionless  Perfect  Gas 

For  an  ideal  gas  with  specific  heats  at  constant  pressure  ( cp )  and  constant  volume 
(cv)  that  do  not  change  with  temperature,  relationships  between  pressure,  temperature, 
density,  Mach  number,  and  other  variables  for  one-dimensional,  frictionless  flow  through 
nozzles  and  diffusers  are  well  developed  (e.g.,  Liepmann  and  Roshko,  1957;  Saad,  1985). 
Those  relationships  ignore  the  weight  of  the  fluid  (i.e.  they  leave  out  the  “p<?”  term  in  Eq. 
(3)  and  the  gdz  term  in  Eq.  (4)).  Because  those  relationships  assume  ideal  gas  behavior, 
they  also  assume  that  no  new  gas  is  being  generated  (for  example,  by  exsolution)  during 
flow.  Dilute  gas/particle  mixtures  in  volcanic  eruptions  have  been  occasionally  modeled 
as  frictionless,  weightless  ideal  gases  (Kieffer,  1981,  1984;  Turcotte  et  al.,  1990).  Such 
models  assume  that  the  erupting  mixtures  roughly  obey  the  ideal  gas  law.  The 
assumption  of  ideal  gas  behavior  tends  to  be  more  valid  as  the  volume  fraction  (or  mass 
fraction)  of  gas  in  the  mixture  increases. 

Using  these  assumptions,  pres  sure- velocity  relationships  of  adiabatically 
decompressing  ideal  pseudogases  follow  the  relationship  (Kieffer,  1984): 


/n,r=constant 


(37) 
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where  y  is  the  ratio  c/cv  of  the  gas/particulate  mixture.  For  air,  7=1.4.  For  FFO  gas,  yis 
generally  lower  (e.g.,  1.236  for  H2O  gas  at  7=900  C,  p= 0.1  MPa).  For  gas/particulate 
mixtures,  the  parameter  yis  calculated  from  the  following  formula: 


7  = 


msCP,g  +  mmCp,m  +  mxCp,x 

mgGv,g  +  +  mxCy,x 


(38) 


By  combining  Eq.  (37)  with  the  continuity  and  momentum  equations  for  an  ideal 
gas,  one  obtains  the  following  relationships  between  pressure  (p ideal ),  density  (p ideal), 
temperature  ( T ideal ),  and  Mach  number  for  flow  within  a  nozzle  (Saad,  1985,  p.  85-88): 


—£—  =  1  +  — — -M2 


ideal 


(39) 


Po 

P  ideal 


7  -1 

1  +  3- - M 

2 


7-1 


(40) 


P, 


ideal 


+i^m=  r 

V  2  V 


(41) 


where  70,  p0,  and  p„  are  the  temperature  (Kelvin),  pressure,  and  density  of  the  mixture  in 
an  upstream  reservoir  where  the  velocity  is  negligible.  If  70,  p0,  and  p0  are  known,  and 
the  Mach  number  at  a  particular  point  in  the  nozzle  is  known,  then  the  temperature, 
pressure,  and  density  at  those  points  can  be  calculated. 

An  ideal  gas/particulate  mixture  can  be  approximated  in  the  program  Conflow  by 
making  the  following  changes:  (1)  set  the  weight  percent  gas  in  the  system  at  100%;  (2) 
set  the  conduit  length  to  be  very  short  to  minimize  the  effects  of  gravity  and  friction  in  the 
calculations.  Flow  through  the  conduit  is  then  calculated  by  setting  a  constant  pressure 
gradient  and  having  the  program  calculate  the  cross-sectional  profile.  The  model 
calculates  the  Mach  number,  temperature,  density,  and  pressure  at  each  point.  Those 
properties  are  plotted  (solid  lines)  as  a  function  of  conduit  position  in  Fig.  9  for  a  5-m 
long  conduit.  At  each  depth,  using  the  Mach  number  calculated  by  Conflow,  the  ideal  gas 
values  of  density,  pressure,  and  temperature  were  calculated  using  Eqs.  (39)-(41).  Those 
values  are  plotted  as  triangles. 

The  ideal  gas  results  are  similar  but  not  identical  to  those  give  by  Conflow. 
Differences  in  the  results  are  assumed  to  be  due  to  (1)  the  non-ideal  properties  of  FFO 
gas,  and  (2)  friction  and  gravity  effects. 
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Installation  and  system  requirements 

The  program  Conflow  can  be  obtained  by  anonymous  ftp  by  pointing  your  web 
browser  to  the  following  USGS  site: 

ftp://elektra.wr.usgs.gov 

Once  entering  this  ftp  site,  go  to  Pub/lgmastin/conflow.  The  program  Conflow  is  in  the 
form  of  a  self-extracting  Zip  file  named  conf  lowzip .  exe.  Source  files  to  the  Fortran 
version  of  this  program  are  in  the  subdirectory  “source  files”.  The  report  you  are  reading 
is  also  available  in  digital  form  at  that  site  as  ofiie.pdf  (in  Portable  Document  Format). 
To  read  the  documentation  file,  you  will  need  Adobe  Acrobat  Reader  ,  which  can  be 
downloaded  free  of  charge  at: 


http://www.adobe.com/products/acrobat/readstep.html 


To  install  the  model,  do  the  following: 

1)  Copy  the  file  "conf  lowzip .  exe"  in  the  above  ftp  directory  to  the  hard  drive  of 
your  Windows-based  computer.  The  zipped  file  occupies  2.2  megabytes  of  disk 
space. 

2)  Double-click  on  the  file  icon.  It  will  then  unzip  and  place  the  unzipped  files  in  a 
new  directory  labeled  "conflow."  The  directory  size  will  be  about  3.1  megabytes. 

3)  Go  into  the  conflow  directory,  and  double-click  on  the  "setup .  exe"  icon.  This 
will  install  the  program,  place  the  executable  files  in  the  directory  "c :  \program 

f  iles\Conf  low",  and  place  the  icon  under  the  "program"  menu  of  the  Start 
button. 

4)  To  Start  the  program,  go  to  Start  »  Programs  »  Conflow. 

If  you  wish  to  uninstall  the  program  later,  you  can  do  so  by  going  to  start  >> 

Settings  >>  Control  Panel,  and  double-clicking  on  the  icon  "Add/ Remove 
Programs".  On  the  "install/uninstall"  tab,  choose  "Conflow"  from  the  list  of 
programs,  then  click  the  "Add /Remove"  command  button. 

Conflow  will  operate  on  any  Windows  -based  computer  running  on  an  Intel  (or 
equivalent)  80386  or  later  processor.  Because  Conflow  is  one-dimensional,  it  does  not 
require  large  amounts  of  memory;  any  recent  Windows  -based  computer  should  be  fast 
enough  to  operate  it.  Informal  tests  using  the  program's  default  input  conditions 
(Kilauean  magma,  1-km  long  conduit)  give  solution  times  ranging  from  about  8  seconds 
on  a  Pentium  II  500  MHz  computer  with  64  MB  RAM  to  about  40  seconds  on  older 
Pentiums  with  about  15  Mb  RAM.  More  silicic  magmas  and  longer  conduits  require 
longer  run  times. 
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Entering  compositional  information 

Launching  the  program  will  open  a  window  (Fig.  10)  which  allows  you  to  choose 
the  composition  of  your  melt,  its  gas  content,  crystal  content,  temperature  and  initial 
pressure.  Here  are  some  tips  on  using  the  window. 
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Figure  1 0:  Magma  composition  page  of  Conflow 

•  Entering  data.  You  can  enter  a  composition  by  typing  in  the  weight  percent  of  the 
constituent  oxides  individually  (in  the  left  column),  or  by  choosing  from  one  of 
several  pre-defined  magma  types  (right  column).  The  program  plots  that 
composition  on  a  silica- alkali  diagram  (right).  The  weight  percent  water,  entered  in 
the  text  box  on  the  lower  left,  refers  to  the  percent  water  in  the  total  mixture 
(melt+crystal+gas),  not  the  dissolved  water  in  the  melt  alone.  The  dissolved  water 
content  is  determined  from  the  total  water  content  and  the  gas  solubility  of  the  melt 
at  the  given  temperature  and  pressure. 

•  Saving  and  loading  data.  You  can  save  compositional  information  to  a  file  by 
choosing  File  >>  save  properties,  or  load  compositional  information  saved 
earlier  by  choosing  File  >>  Load  properties."  Files  of  compositional  data  are 
in  ASCH  format  and  have  the  suffix  “.mpr”. 
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Figure  1 1 :  windows  for  thermodynamic  properties,  water  solubility,  and  viscosity  of  melt-gas- 
crystal  mixture 

•  Choosing phenocryst  types.  In  the  lower  box,  you  can  also  choose  the  dominant 
type  of  phenocryst  in  the  melt  and  its  abundance  (as  a  volume  percent  of  the  liquid- 
crystal  mixture).  By  choosing  a  phenocryst  type  using  the  radio  buttons,  the  program 
will  calculate  its  specific  heat  and  density  using  relations  from  Berman  (1988). 

•  You  can  view  the  thermodynamic  properties ,  the  gas  solubility  in  the  melt,  and  the 
mixture  viscosity  by  clicking  the  appropriate  command  buttons.  These  command 
buttons  will  open  additional  windows  (Fig.  1 1)  with  plots  and  output  information. 
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Figure  1 2:  Window  for  conduit  properties 

Specifying  conduit  properties 

By  clicking  on  the  "conduit  properties"  button  in  the  Magma  composition  window, 
a  new  window  opens  (Fig.  12)  from  which  you  can  specify  properties  of  the  conduit.  The 
following  is  a  list  of  variables  and  their  effects  on  the  model.  Their  effects  are  described 
in  greater  detail  in  the  section  "RUNNING  THE  MODEL  FROM  THE  COMMAND 
LINE". 

Execution  options.  The  conduit  model  calculates  flow  properties  in  the  conduit 
using  one  of  two  assumptions:  (1)  The  cross-sectional  area  of  the  conduit  is  specified  as 
input  to  the  program,  and  the  program  determines  the  pressure  profile;  or  (2)  the  pressure 
profile  is  specified  as  input,  and  the  model  finds  the  conduit  geometry  that  produces  that 
pressure  profile.  By  clicking  on  one  of  the  radio  buttons  in  the  Conduit  Properties  option 
box  on  the  upper  right,  you  are  choosing  among  those  options. 

Iteration  control.  Normally,  if  the  conduit's  cross  sectional  area  is  specified,  the 
conduit  adjusts  the  input  velocity  at  the  base  of  the  conduit  until  either  (1)  the  velocity  at 
the  top  of  the  conduit  equals  the  sonic  velocity  (choked  flow),  or  (2)  the  pressure  at  the 
top  of  the  conduit  equals  that  specified  in  the  text  box  of  this  window.  By  checking  the 
"ignore  upper  boundary  conditions"  radio  button  in  the  Iteration  Control  option  box,  the 
program  will  ignore  the  upper  boundary  conditions;  it  will  simply  calculate  a  single  run 
up  the  conduit,  using  the  input  velocity  provided. 
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Temperature  control.  Conflow  is  capable  of  calculating  adiabatic  temperature 
changes  within  the  conduit  as  a  result  of  shear  heating,  gas  expansion,  and  gas  exsolution. 
By  clicking  the  "Constant  temperature"  radio  button  in  the  Temperature  Control  option 
box,  you  can  convert  the  program  to  isothermal  calculations.  Isothermal  calculations  are 
somewhat  faster  than  those  that  consider  adiabatic  temperature  changes. 

Depth,  pressure,  and  gravitational  constant.  Normally,  the  pressure  and  depth  at 
the  top  of  the  conduit  are  set  to  0.1013  MPa  (1  atm)  and  0  meters,  respectively;  but  you 
can  change  these  if  you  prefer  to  model  only  a  section  of  the  conduit,  not  ending  at  the 
ground  surface.  Similarly,  the  gravitational  constant  and  pressure  at  the  top  of  the 
conduit  can  be  changed  to  model  eruptions  on  other  planets.  The  depth  at  the  base  of  the 
conduit  can  be  considered  the  depth  immediately  above  a  magma  chamber,  though  this 
assumption  is  not  necessary.  Any  depth  above  a  magma  chamber  can  be  used  as  a 
starting  point  for  the  model.  The  pressure  at  the  base  of  the  conduit  is  considered  by 
many  modelers  to  be  near  the  lithostatic  pressure  at  that  depth  (the  lithostatic  pressure 
gradient  is  usually  about  20-25  MPa/km).  In  reality,  the  magma  pressure  could  vary  by 
tens  of  percent  from  the  lithostatic  pressure  depending  (among  other  things)  on  the  degree 
of  anisotropy  of  in  the  principle  stresses  of  the  host  rock,  the  conduit  shape,  and  the  rock 
strength. 

Input  velocity.  In  cases  where  the  user  specifies  the  conduit  geometry  and  requires 
that  the  upper  boundary  condition  be  satisfied,  the  program  uses  the  input  velocity 
specified  in  this  text  box  as  the  starting  point  of  an  iterative  sequence.  In  successive 
model  runs  it  adjusts  the  input  velocity  until  either  (1)  the  velocity  at  the  top  of  the 
conduit  is  sonic  (choked  flow);  or  (2)  the  pressure  at  the  top  of  the  conduit  matches  the 
final  pressure  specified. 

Name  of  output  file .  The  program  will  generate  a  long  output  file  whose  name  can 
be  specified  in  this  text  box.  The  file  will  be  in  normal  ASCII  format  but  needn't  have  the 
.txt  suffix  that  designates  it  as  a  text  file. 

Variables  for  output.  By  clicking  on  this  command  button,  you  open  a  new 
window  in  which  you  can  specify  the  calculated  flow  properties  (up  to  seven)  to  be 
written  to  the  output  file,  and  which  properties  (up  to  four)  will  be  plotted.  If  you  have 
already  run  previous  models  during  this  programming  session,  this  button  will  be 
disabled  so  that  output  variables  will  be  consistent  from  one  run  to  another. 

Running  the  model. 

By  clicking  the  "Run  model"  command  button,  you  run  the  numerical  model  for 
conduit  flow,  using  the  input  values  that  were  defined  in  this  window  and  the  Magma 
composition  window.  The  numerical  model  opens  a  DOS  window  (Fig.  13)  and  writes 
out  intermediate  results  to  the  screen  as  it  determines  a  solution.  The  meaning  of  the 
information  written  to  the  screen  during  execution  is  described  in  the  section  "Model 
Execution." 


Using  the  Windows-based  Version 


27 


START [N«  him  won up  9  m 

ruo  t  f  lnw 

-  H.lUMlF+i 

Ah  k%r'Z 

i 

p  taJ 

ft  cirlPdS 

r-id  j  ns 

Ut  H'.lft  Li:  1 

trtrj 

uian 

anrli  R 

i 

-i  HGfl.HiiH 

2  h  .000 

5  .  00S1 

a'm? 

i  .am 

1  .959 

0  .  tlHt 

m 

-H.2F1 

1 .47S 

S.ffiW 

a.  ns  ? 

]  .051 

1  .938 

0.001 

m 

-62.6fl» 

]  .325 

5  .  OflH 

a.2»H 

1  .262 

1.757 

0.001 

Ml 

-EE.  610 

1 . 170 

5 .  m 

0.333 

]  .500 

1.699 

H.  002 

n 

47. 992 

1  .067 

5  .mm 

0."]3H 

1  .776 

1 .631 

0 .  002 

1«« 

-39.835 

0.955 

e  .0«a 

0.525 

2.103 

1.553 

0.002 

12H 

31.  M? 

0.854 

5 . 11I.1VI 

0.575 

2.491 

1.4V5 

0 . 003 

140 

21.737 

0.763 

5 . 0«a 

0.55V 

2.727 

1.427 

0 . 001 

L6tl 

11 .920 

0.680 

5.00a 

0.711 

3.446 

1 .162 

0 . 004 

19*1 

-1 .31* 

H.60E 

s  .n«n 

0.754 

4.K53 

1 .266 

0 . 005 

193 

Ol.titJH 

0.577 

5 .0tWI 

0.75V 

4.133 

1.269 

0.005 

exit  pniiure  > 

pF  In  al  aim] 

H  C  1 

try  lit  g 

It  it', p  inyiil 

vr  in  city 

t.  bmsau 

fi-i'e 

STAJtttHC  HUH  HIWKO  2- 

htd.BE  f  1UK 

(1 .31 641:  +0Ft  Jtsj/E 

i 

a  <m> 

p  (mp*? 

rdt]  luf 

(if  IMS  (IB  1 

<WS> 

1  it  y 

(fl®  5 

Hiielt  h 

1 

-10W.UBH 

Z6  . uwn 

5.HIW 

0  (10(1 

1  .500 

1  K59 

0.002 

"69,714  4 

zf r  8  mTTC  VvtB 

5  ,  (Uftl 

a.  055 

i.bnv 

1  ,s:i9 

U.083 

4M 

-ft  3-  rj.  7 

1  .32? 

5  .  HUU 

H, 

Figure  13:  DOS  window  that  opens  when  the  conduit  model  is  launched. 
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Figure  14:  Window  displaying  output  to  model. 
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Viewing  output 

Once  the  model  has  finished  and  the  DOS  window  has  closed,  the  output  will  be 
written  as  text  to  the  file  given  in  the  "name  of  output  file"  text  box.  You  can  view  the 
output  data  in  tabular  form  by  clicking  the  "view  output"  command  button.  It  will  bring 
up  a  new  window  (Fig.  14)  showing  a  summary  of  the  input  conditions  and  results  of 
iterative  model  runs  shown  in  the  upper  text  box,  and  a  table  of  the  final  flow  properties, 
as  a  function  of  depth,  in  the  lower  spreadsheet. 

Plotting  output 

By  pressing  the  plot  command  button  in  the  results  window,  a  new  window 
(Fig.  15)  will  appear  with  plots  of  the  variables  that  were  chosen  in  the  Output 
variables  window,  accessed  through  the  command  button  on  the  input  Properties 
window.  If  you  executed  any  other  model  runs  since  opening  Conflow,  those  model  runs 
will  also  be  plotted  for  comparison.  An  explanation  of  the  abbreviations  used  for  x-axis 
labels  is  provided  under  Help  >>  more  info.  You  can  label  the  plot  and  print  it  out  if 
you  wish. 
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Figure  15:  Window  displaying  plotted  output. 

Running  the  M  odel  from  the  Command  Line 

The  graphical  and  interactive  version  of  this  program  that  runs  on  Windows-based 
computers  calls  a  simple  Fortran  program  (named  conf  ort .  exe)  from  the  DOS 
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command  line.  The  program  confort.exe  resides  in  the  same  directory  as  Conflow.exe 
(usually  c:\program  files\Conflow).  Confort.exe  can  also  be  operated  from  a  DOS 
command  line  by  opening  a  DOS  window,  moving  to  the  directory  containing 
confort.exe,  and  typing  "confort".  Compiled  versions  of  this  program  will  also  be  posted 
in  the  ftp  directory  Elektra.wr.usg.gov/Ftp_Access/Pub/lgmastin/Conflow,  which  will  run 
on  other  operating  systems. 

When  executed,  the  program  reads  from  the  ASCII  input  file  conin,  which  resides 
in  the  same  directory  and  can  be  edited  using  any  text  editor.  The  file  appears  as  follows: 


INPUT  PARAMETERS: 
hawaii . out 
diam 

27.,  0.1013 
2 

1.0 

1145.,  2 
1000  . 

0.27 

2 

1000  ,  0. 

5.  ,  5. 

9.81 

51 . 96 
14.21 
0.00 

10.96 
6.59 
10.86 
2.53 
2.48 

0 .416 
0.  2600 


PARAMETER  EXPLANATIONS: 
name  of  output  file 

specify  conduit  diameter  (diam)  or  pressure  gradient  (pgrd) 
initial,  final  pressure  (MPa) 
iteration  number* 
initial  velocity  (m/s) 

initial  temperature  (C) ,  itemp  (l=const  T,  2=variable  T) 

specific  heat  of  crystals  (J/kg  K) 

h2o  content  (wt%) 

vesiculation  parameter** 

initial,  final  depth  (m) 

conduit  diameter  (m)  at  bottom,  at  top 

gravitational  acceleration  (m/s2) 

wt%  Si02  (anhydrous) 

wt%  A1203  (anhydrous) 

wt%  Fe203 

wt%  FeO 

wt  %  MgO 

wt%  CaO 

wt  %  Ti02 

wt  %  Na20 

wt  %  K20 

volume  %  crystals,  xtl  density  (kg/m3) 


NOTES  ON  INPUT  PARAMETERS: 

‘iteration  #=2  if  the  velocity  is  to  be  adjusted  automatically  to  reach 
sonic  velocities  at  the  exit  (valid  only  if  icalc=l),  or 

1  if  no  adjustment  is  desired. 

“vesiculation  p.=  2  if  gas  exsolution  is  to  stop  after  fragmentation 
1  if  not 


Output  Parameters: 

List  of  variables  to  be  written  out.  Enter  a  number  in  the 

first  column  indicating  the  column  #  where  this  variable  will  be  written 

in  the  output  file.  You  can  write  out  up  to  seven  variables, 
x-sectional  area  (m2) 

7  Mach  number 

5  pressure  (MPa) 

3  log  Reynolds  number 
mixture  density 

time  (s)  since  entering  conduit 

6  velocity  (m/s) 

4  volume  fraction  gas 
log  viscosity  (Pa  s) 

1  z  (depth,  meters) 
d(x-s  area)/dz,  meters 
log  pressure  (MPa) 

dpdz  (pressure  gradient.  Pa/m) 
log  dz  (vert,  step  size,  m) 
f  (friction  factor) 
gamma  (Cp/Cv  for  gas  phase) 
mf  (mass  fraction  exsolved  gas) 
mm  (mass  fraction  magma) 
r  (Universal  Gas  const.  *  n) 
rhof  (gas  density) 
sv  (sonic  velocity  (m/s) 
temperature  (C) 
enthalpy  of  mixture  (kj/kg) 
cp  (sp.  heat)  of  gas  (kj/kg  C) 

2  conduit  radius (m) 
dissolved  h2o  (wt%) 
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cp  (sp.  heat)  of  melt  (kj/kg  C) 


The  twenty-one  lines  following  the  first  line  of  the  file  contain  the  input  parameters  on 
the  left  side.  Those  parameters  are  read  using  unformatted  read  statements,  so  they  can  be 
changed  without  worrying  about  column  numbers  or  number  of  decimal  places.  Just  be 
careful  not  to  add  or  delete  any  lines  while  editing  the  file.  All  variables  are  double 
precision,  real  numbers,  with  the  exceptions  of  the  vesiculation  parameter,  the  iteration 
number,  and  itemp,  which  are  integers,  and  the  parameter  on  the  second  line,  which  is  a 
4-character  variable  ("diam"  or  "pgrd"). 

The  right-  hand  side  of  each  line  explains  (briefly)  what  each  parameter  represents. 
Parameter  explanations  that  require  somewhat  more  information  are  followed  by 
asterisks,  with  supplemental  information  on  following  lines.  Although  most  parameters 
are  self-explanatory,  the  following  parameters  require  more  detailed  information: 

Specifying  conduitdiameteror  pressure  gradient 

The  second  line  of  the  input  file  specifies  which  option  to  use  when  running  the 
program.  If  "diam"  is  specified,  the  program  assumes  a  constant  conduit  diameter  and 
calculates  a  pressure  profile.  If  "pgrd"  is  specified,  constant  pressure  gradient  is  assumed 
and  the  program  calculates  the  profile  in  cross-sectional  area  that  would  produce  such  a 
pressure  gradient. 

Pressure  at  base  and  top  of  conduit 

This  parameter  is  used  only  if  the  conduit  diameter  is  specified  (instead  of  the 
pressure  gradient).  There  is  no  real  upper  limit  to  the  maximum  input  pressure  that  can  be 
used,  but  the  lower  limit  is  constrained  by  the  weight  (per  unit  area)  of  the  magma  in  the 
conduit.  If  the  input  pressure  is  less  than  that  weight  per  unit  area,  the  magma  will  not 
erupt.  In  such  a  case,  the  model  will  reach  p=l  atm  at  some  depth  below  the  surface.  If 
the  model  is  set  to  iterate  until  p=l  atm  or  M=  1  at  the  surface,  it  will  decrease  the  velocity 
at  the  base  of  the  conduit  and  try  another  run.  If,  after  several  iterations,  the  initial 
velocity  drops  below  0.001  m/s  and  p=l  atm  is  still  reached  below  the  ground  surface,  the 
program  returns  the  following  message  to  the  screen: 

pressure  insufficient  to  produce  eruption 

and  writes  the  results  of  the  last  run  (in  which  initial  velocity=0.001  m/s)  to  the  output 
file.  The  following  table  indicates  the  minimum  pressures  that  will  produce  upflow  for 
various  conduit  lengths  for  Kilauean  basalt,  given  other  input  parameters  shown  in  the 
example  input  file: 


depth  at  base  of  conduit 

minimum  pressure  for 

100  m 

0.17 

200 

0.35 

500 

4.35 

1000 

17.1 

3000 

71.6 

Running  the  M  odel  from  the  Command  Line 


31 


For  silicic  melts,  minimum  pressures  that  will  still  produce  an  eruption  are  negligible 
(less  than  5  MPa  at  the  base  of  a  5-km  long  conduit  for  Pinatubo  magma,  for  example; 
Fig.  12).  In  reality,  significantly  higher  pressures  would  be  necessary  at  these  depths  to 
drive  eruptions,  since  gas  escape  at  low  magma  velocities  would  densify  the  magma 
column  and  increase  its  weight. 

If  the  combination  of  input  pressure  and  FFO  content  of  the  melt  are  such  that 
vesicularity  at  the  base  of  the  conduit  exceeds  75%,  the  program  assumes  that  the  melt 
has  already  fragmented  before  entering  the  base  of  the  conduit. 

In  general,  the  pressure  at  the  top  of  the  conduit  is  specified  to  be  1  atmosphere 
(0.1013  MPa).  However  this  program  gives  the  user  the  option  to  specify  other  output 
pressures,  which  may  be  useful  under  three  circumstances:  (1)  when  modeling  eruption 
dynamics  on  other  planets;  (2)  when  modeling  magma  flow  through  a  particular  depth 
interval  whose  top  does  not  lie  at  the  surface;  or  (3)  modeling  an  eruption  that  vents  to  the 
sea  floor,  the  floor  of  a  lake,  or  of  a  lava  pond.  The  option  of  variable  final  pressure 
allows  the  user  to  model  a  complicated  conduit  geometry  by  breaking  the  conduit  into 
sections  and  modeling  each  depth  interval  separately. 

If  option  2  (pressure  gradient  specified)  is  chosen,  the  pressure  within  the  conduit  is 
assumed  to  vary  linearly  between  the  initial  and  final  values.  The  pressure  gradient  is  the 
difference  between  these  pressures,  divided  by  the  length  of  the  conduit.  The  conduit's 
cross-sectional  area  is  adjusted,  along  with  flow  properties,  to  fit  this  gradient.  Some 
models  of  conduit  flow  (e.g.,  Wilson  and  Head,  1981;  Dobran,  1992)  assume  that  the 
pressure  gradient  driving  magma  flow  is  the  gradient  paQ,  determined  by  the  country 
rock  density,  pcr.  In  those  programs,  if  a  country  rock  density  of  2300  kg/m3  is  used  as 
input,  the  program  calculates  a  pressure  gradient  of  pcrg=  2.25xl04  Pa/m,  and  a  pressure 
at  the  base  of  a  3-km-long  conduit  of  1.013xl05  Pa  +  (3000m)(2.25xl04  Pa/m)  = 

6.78x10  Pa,  or  67.8  MPa.  In  fact,  far-field  horizontal  stress  gradients  may  be  as 
important  as  the  lithostatic  pressure  gradient  in  controlling  the  flow  up  the  conduit.  In  the 
program  Conflow,  pressures  at  the  top  and  bottom  of  the  conduit  are  given  directly  as 
input  to  the  program  rather  than  a  rock  density  from  which  a  pressure  gradient  is 
calculated. 

There  is  one  caveat  when  considering  the  input  value  for  pressure.  If  the  pressure 
gradient  in  the  conduit  is  less  than  that  due  to  the  weight  of  the  magma-gas  mixture  at  the 
base  of  the  conduit,  the  magma  may  not  flow  upward.  In  that  case,  the  following  error 
message  will  appear: 

WARNING! ! 

Density  of  magma/gas  mixture  =  2294.3  kg/m3 . 

Thus  its  pressure  gradient  is  22.9  MPa/km. 

This  is  greater  than  that  specified  for  the  conduit:  22.9 

The  conduit  diameter  will  probably  grow  to  an 
unrealistically  large  value  before  reaching 
the  surface.  Do  you  wish  to  continue  (y/n) ? 


The  program  will  increase  the  in  cross-sectional  area  with  depth  to  match  the  specified 
pressure  gradient;  but  the  pressure  gradient  nevertheless  not  be  reached  without 
expanding  the  conduit  diameter  to  an  unrealistically  large  value.  If  this  happens,  you  will 
see  the  following  error  message: 
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Conduit  diameter  or  dadz  is  unrealistically  large. 

Program  stopped. 

You  will  have  to  increase  the  pressure  gradient  (by  either  increasing  the  initial  pressure  or 
decreasing  the  final  pressure)  and  try  again. 

Iteration  number 

If  this  number  is  1  and  the  conduit  diameter  (rather  than  pressure  gradient)  is 
specified,  the  velocity  will  not  be  adjusted  to  match  the  exit  boundary  conditions.  The 
program  calculates  a  single  run  up  the  conduit  and  writes  out  the  results  without 
attempting  to  match  the  exit  conditions  to  the  boundary  conditions.  If  the  velocity  of  the 
mixture  reaches  the  sonic  velocity  before  the  calculations  reach  the  top  of  the  conduit,  the 
program  stops  at  that  depth.  The  same  is  the  case  if  the  pressure  drops  below  atmospheric 
before  the  calculations  reach  the  surface. 

Tolerance  levels.  If  the  iteration  number=2,  the  program  will  iterate  until  the 
output  pressure  is  between  0.1012  and  0.1014  MPa  (1  atm=  0.1013  MPa).  For  the  M=  1 
boundary  condition,  the  program  iterates  until  M=1  is  reached  (to  double -precision 
accuracy)  less  than  2  meters  below  the  surface. 

Lack  of  Convergence.  On  a  few  occasions,  the  program  may  have  some  difficulty 
reaching  a  solution  within  the  tolerance  levels  specified  above.  Sometimes  this  problem 
is  due  to  the  fact  that  final  exit  pressures  or  velocities  are  extremely  sensitive  to  the  input 
velocity,  and  very  slight  changes  in  input  velocity  (usually  less  than  10  4  m/s)  cannot 
produce  an  acceptable  result.  In  such  a  case,  the  program  stops,  writes  out  the  results  of 
its  best  run,  and  prints  the  following  message  to  the  screen: 


limit  of  resolution  reached 

On  more  rare  occasions,  the  program  just  won’t  converge  at  all.  If  this  happens,  a  slight 
change  to  an  input  parameter  will  usually  solve  the  problem. 

Initial  velocity 

In  option  1  (where  conduit  diameter  is  specified),  if  the  iteration  numbcr=l ,  the 
velocity  is  adjusted  until  the  output  pressure=l  atm  or  the  output  velocity=sonic  velocity 
of  the  mixture.  Under  these  circumstances,  the  initial  input  velocity  is  only  the  starting 
point  of  the  iteration  sequence.  If  option  2  is  specified,  or  the  iteration  number=2,  then 
the  initial  velocity  is  used  for  the  final  solution. 

Initial  temperature 

Used  to  calculate  viscosity  of  magma,  specific  volume  of  the  gas  phase  (using  the 
Haar  et  al.  gas  relationships),  and  enthalpy  of  the  magma-gas  mixture.  The  initial 
temperature  is  given  in  degrees  Celsius. 

H2O  content 

This  is  the  amount  of  H2O  (both  dissolved  and  exsolved)  in  the  erupting  mixture,  in 
weight  percent,  we  have  successfully  used  water  contents  from  0%  to  100%,  though  the 
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model  is  not  designed  to  handle  those  end  members.  Real  values  of  gas  entrained  in  melt 
during  plinian,  sub-plinian,  or  lava-fountain  emptions  rarely  exceed  several  percent. 

Vesiculation  parameter 

In  the  uppermost  several  tens  of  meters  of  the  conduit,  when  vesicularity  and 
eruptive  velocities  are  high,  gas  exsolution  rates  may  not  keep  pace  with  the  rate  of 
depressurization.  If  the  vesiculation  parameter  is  set  to  2,  gas  exsolution  is  no  longer 
computed  once  the  magma  fragments  (though  gas  expansion  due  to  decompression  is  still 
calculated).  If  the  vesiculation  parameter  1,  gas  exsolution  is  also  calculated  at  all 
vesicularities.  At  the  base  of  the  conduit,  the  dissolved  water  content  of  the  melt  is 
assumed  to  be  at  equilibrium,  whether  the  magma  has  fragmented  or  not. 

Initial,  final  depth 

The  depth  of  the  base  and  top  of  the  conduit,  in  meters.  Numerically,  the  program 
can  handle  any  arbitrary  starting  depth,  from  several  kilometers  (or  more),  essentially  up 
to  the  ground  surface.  If  unusually  shallow  starting  depths  are  used,  the  mixture  will 
already  be  highly  vesiculated,  possibly  fragmented.  If  the  vesicularity  at  the  base  of  the 
conduit  exceeds  75%,  the  program  assumes  that  the  melt  has  already  fragmented. 

Conduitdiameter  at  base,  attop 

Under  option  1  (where  conduit  diameter  is  specified),  these  variables  give  the  diameter  of 
the  conduit  at  the  base  and  top,  respectively.  The  program  assumes  that  the  conduit's 
diameter  varies  linearly  between  these  values.  If  option  2  (pressure  gradient)  is  specified, 
the  program  reads  only  the  diameter  at  the  base  of  the  conduit.  The  diameter  at  the  top  is 
calculated. 

Gravitational  acceleration 

This  is  normally  9.81  m/s  for  eruptions  on  Earth.  It  can  be  adjusted  for  eruptions 
on  other  planets. 

M  agma  composition 

The  melt  composition  is  given  in  weight  percent  of  the  major  oxides  in  an 
anhydrous  melt.  These  values  are  used  to  calculate  viscosity,  density,  gas  solubility, 
specific  heat,  and  enthalpy  of  the  melt.  The  relations  that  calculate  those  properties  are 
calibrated  for  the  range  of  silicate  melts  found  in  nature.  It  is  not  known  how  valid  these 
properties  are  when  extended  outside  the  range  of  natural  silicate  melts.  When  reading 
the  compositional  data,  the  program  checks  to  see  that  they  add  up  to  100%.  If  the  total 
differs  by  more  than  0.1%  from  100%,  you  will  receive  the  following  error  message: 


total  of  component  oxides  does  not  equal  100%. 
Adjust  automatically?  (y/n) 


If  you  enter  'y',  the  program  will  adjust  the  weight  percent  of  each  oxide  proportionately 
so  that  the  total  equals  100%.  If  you  enter  'n',  the  program  will  stop  and  you  must  edit  the 
input  file  and  try  again. 
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Specifying  the  Variables  to  be  Written  as  Output 

The  last  37  lines  of  the  input  file  contain  the  names  of  variables  that  can  be  written 
to  the  output  file  for  each  depth.  You  must  specify  seven  variables  to  be  written  out.  For 
each  variable  to  be  written  out,  enter  a  number  at  the  beginning  of  the  appropriate  line 
corresponding  to  the  column  in  the  output  file  where  this  variable  will  appear.  In  the 
example  input  file  above,  the  depth  (z)  is  to  be  written  in  the  first  column  of  the  output 
table;  temperature  in  the  second,  velocity  in  the  third,  log  viscosity  in  the  fourth,  and  so 
on. 

You  must  specify  output  variables  for  seven  columns  of  output.  If  you  neglect  to 
specify  output  for  a  given  column  (e.g.,  column  1),  you  will  receive  an  error  message  like 
the  following: 


You  have  entered  your  output  variables  incorrectly.  They  are: 
column  1  not  given 
column  2  temp(C) 
column  3  vel  (m/s) 
column  4  log  vise 
column  5  vfgas 

column  6  p  (MPa) 
column  7  mach  # 

Please  edit  the  input  file,  entering  seven  output  variables,  and  start  again. 
Remember  to  enter  the  output  column  number  in  the  FIRST  column  of  each  appropriate 
line . 


You  can  also  receive  this  error  message  if  the  column  number  specified  is  not 
entered  as  the  FIRST  character  in  the  appropriate  line  of  the  input  file  (i.e.  if  the  column 
number  is  preceded  by  one  or  more  spaces).  Alternatively,  if  you  specify  the  same 
column  number  for  two  or  more  output  variables,  you  may  receive  an  error  message  as 
follows: 


Output  to  column  1  has  been  specified  for  TWO  variables: 
z  (m) 
temp (C) 

Please  edit  the  input  file  to  correct  the  problem. 
Remember  to  enter  the  output  column  number  in  the 
FIRST  column  of  the  appropriate  line  in  the  input  file. 


M  odel  Execution 

The  program  can  be  executed  by  moving  to  the  directory  where  it  resides  and  typing 
“confort”  on  the  command  line.  (If  your  computer  uses  Microsoft  Windows  ,  you 
should  open  a  DOS  window  before  executing  the  program).  As  long  as  the  input  file 
conin,  is  in  the  same  directory  as  the  executable  file,  the  program  should  be  able  to  find 
it.  Two  examples  of  program  execution  are  given  below:  one  using  option  1  (specified 
conduit  diameter),  the  other  using  option  2  (constant  pressure  gradient). 

Example  using  option  1 

Once  the  program  is  started,  it  will  write  out  the  input  values  to  the  DOS  window  as 
follows: 


using  program  "Conflow" 


5  The  use  of  trade  names  is  not  intended  to  be  an  endorsement  of  those  products. 
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INPUT  VALUES: 


input  velocity  =  1.0000  m/s 

magma  density=  2666.  kg/m3  (calculated) 
input  temperature  =  1145.  degrees  Celsius 
initial  dissolved  h2o=  0.270  wt  % 
depth  at  base  of  conduit=  1000.0000  m 
depth  at  top  of  conduit=  0.0000  m 

MAGMA  COMPOSITION: 


51 . 960 

SI02  (wt%  anhydrous) 

14.210 

A1203 

0.000 

Fe203 

10.960 

FeO 

6.590 

MgO 

10.860 

CaO 

2.530 

Ti02 

2.480 

Na20 

0.416 

K20 

0.000 

crystal  content  (vol%  of  melt) 

2600.000  crystal  density  (kg/m3) 


calculating  temperature  change 

specified  conduit  diameter: 
diameter  at  base  =  5.00  meters 

diameter  at  top  =  5.00  meters 

input  pressure  =  27.00  MPa 

automatic  velocity  adjustment 

no  exsolution  after  fragmentation 

vfgas=0.75  is  fragmentation  criterion 


These  are  the  same  input  parameters  specified  in  the  sample  input  file  above.  For  this 
run,  the  conduit  diameter  is  taken  to  be  constant  and  the  program  adjusts  the  input 
velocity  until  M=  1  or  p=  1  atm  at  the  surface. 

Next,  the  program  will  begin  calculating  flow  properties  from  the  bottom  to  the  top 
of  the  conduit.  The  output  to  the  screen  at  this  point  in  execution  is: 


STARTING  RUN  NUMBER  1: 


mass  flux=  0.5244E+05  kg/s 


i  z  (m)  radius  log  Re 

1  -1000.000  2.500  2.319 

92  0.000  2.500  2.014 

exit  pressure  >  pfinal  and  M  <  1 


vfgas  p  (MPa)  vel  (m/s)  mach  # 
0.000  27.000  1.000  0.000 

0.504  0.936  2.014  0.054 


After  writing  out  the  mass  flux  calculation,  the  program  has  written  a  line  of  output 
variables  calculated  at  the  bottom  of  the  conduit,  and  a  second  line  at  the  final  depth.  On 
the  last  line,  the  program  notes  that  the  final  exit  pressure  exceeds  1  atm  and  the  Mach 
number  is  less  than  1.  The  program  therefore  increases  the  input  velocity  and  computes  a 
second  run,  writing  the  output  as  follows: 


trying  new  input  velocity  1.50000  m/s 


STARTING  RUN  NUMBER  2: 


mass  flux=  0.7867E+05  kg/s 


i  z  (m)  radius  log  Re 

1  -1000.000  2.500  2.495 

95  0.000  2.500  2.175 

exit  pressure  >  pfinal  and  M  <  1 


vfgas  p  (MPa)  vel  (m/s)  mach  # 
0.000  27.000  1.500  0.000 

0.524  0.912  3.149  0.085 
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Again,  the  final  exit  velocity  exceeded  1  atm  and  the  final  Mach  number  is  less  than  1. 
After  adjusting  the  input  velocity  again,  a  third  run  is  attempted: 


trying  new  input  velocity 

STARTING  RUN  NUMBER  3: 

i  z  (m)  radius 

1  -1000.000  2.500 

100  0.000  2.500 

exit  pressure  >  pfinal  an 


2.25000  m/s 

mass  flux= 

log  Re  vfgas  p 

2.671  0.000 

2.329  0.555 

M  <  1 


0 . 1180E+06  kg/s 

(MPa)  vel  (m/s)  mach  # 
27.000  2.250  0.000 

0.873  5.050  0.139 


After  increasing  the  input  velocity  three  more  times,  Conflow  exceeds  Mach  1  before 
reaching  the  surface: 


STARTING  RUN  NUMBER  6: 


mass  flux=  0.3982E+06  kg/s 


i  z  (m) 

1  -1000.000 
244  -12.930 

mach  number  > 


radius  log  Re  vfgas  p  (MPa)  vel 

2.500  3.200  0.000  27.000 

2.500  8.985  0.773  0.523 

1.  adjusting  initial  velocity 


(m/s) 

7.594 

33.334 


mach  # 
0.001 
1.000 


After  several  more  adjustments  to  the  input  velocity,  Conflow  final  reaches  an  acceptable 
solution: 


trying  new  input  velocity  6.50391  m/s 


STARTING  RUN  NUMBER 


mass  flux=  0.3411E+06  kg/s 


i  z  (m)  radius  log  Re 

1  -1000.000  2.500  3.132 

257  -0.082  2.500  8.959 


vfgas  p  (MPa)  vel  (m/s)  mach  # 
0.000  27.000  6.504  0.001 

0.799  0.447  32.234  1.000 


successful  completion 


AFTER  ISENTROPIC  EQUILIBRATION  TO  1  ATM  PRESSURE: 
final  temperature  =  1143.19  deg.  C 
temperature  change  =  0.691  deg.  K 

enthalpy  change  =  0.1117E+04  J/kg 

max.  theoretical  velocity  =  79.49  m/s 


maximum  water  table  depth  that  will  allow  g.w.  influx  =  34.53  meters 

Negative  values  are  below  the  ground  surface, 
positive  values  are  above. 


This  output  shows  that,  during  the  last  run,  the  Mach  number  reached  1  at  a  depth  of 
0.082  m,  well  within  the  tolerance  limit  of  2  m. 

In  all  runs  where  the  Mach  number=l  when  the  mixture  exits  the  conduit,  the 
pressure  will  be  greater  than  1  atmosphere.  After  the  mixture  leaves  the  conduit,  it  will 
continue  to  accelerate  and  cool  adiabatically  as  it  drops  to  atmospheric  pressure.  If  we 
assume  that  these  processes  take  place  isentropically  (i.e.  without  friction),  we  can 
calculate  a  maximum  theoretical  velocity  and  a  maximum  amount  of  adiabatic  cooling. 
These  calculations  are  done  by  assuming  that  all  excess  enthalpy  in  the  mixture  is 
converted  to  kinetic  energy  during  expansion  (Mastin,  1995a).  Procedures  for  this 
calculation  are  explained  in  Appendix  D.  The  output  written  above  indicates  that  the 
velocity  could  theoretically  accelerate  from  32.234  m/s  to  79.49  m/s  after  leaving  the 
vent.  During  expansion,  the  mixture  would  theoretically  cool  by  about  a  half  degree 
Celsius. 
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A  final  calculation  is  made  of  the  depth  of  the  water  table  required  to  produce 
groundwater  influx  during  the  eruption.  The  calculation  is  made  by  numerically  drawing 
a  hydrostatic  pressure  curve  that  is  just  tangent  to  the  pressure  profile  in  the  conduit.  The 
depth  (“WT”  in  Fig.  16)  at  which  the  hydrostat  reaches  one  atmosphere  gives  the  water 
table  depth  listed  above,  //'subsurface  water  pressures  follow  the  hydrostatic  curve,  then 
a  water  table  at  this  depth  or  higher  would  create  hydrostatic  pressures  sufficient  to  drive 
water  into  the  conduit.  Whether  water  enters  in  sufficient  quantities  to  produce 
explosive,  phreatomagmatic  interactions,  depends  on  factors  such  as  rock  permeability, 
that  are  not  considered  here.  This  model  also  does  not  consider  other  important  processes 
that  take  place  once  water  enters  a  conduit  and  mixes  with  magma,  possibly  including 
steam  expansion,  brittle  fragmentation  of  melt  caused  by  high  strain  rates,  and  fuel- 
coolant  interactions. 


Figure  1 6:  flow  properties  in  a  5-km  long,  50-m  diameter  eruptive  conduit  containing  Pinatubo 
magma  at  780°  C,  with  6  wt%  water,  for  three  different  input  pressures  at  the  base  of  the  conduit. 
The  plot  illustrates  (a)  that  extremely  low  input  pressures  can  still  generate  eruptions,  and  (2)  that 
the  final  eruptive  velocity  is  extremely  insensitive  to  changes  in  input  pressure.  The  heavy  solid 
line  illustrates  a  hydrostatic  gradient  that  is  just  tangent  to  the  solid  pressure  profile.  The  y  value 
of  this  line  at  p=1  atm  (labeled  “WT”)  gives  the  maximum  water-table  depth  of  a  normally 
pressured  groundwater  system  in  which  water  pressures  would  be  sufficient  to  drive  water  into  the 
conduit. 

Program  output.  Once  the  program  completes  its  calculations,  open  the  output  file 
and  you  will  see  the  following  table  (already  described): 

OUTPUT  TABLE 

257  i  z  (m)  radius  log  Re  vfgas  p  (MPa)  vel  (m/s)  mach  # 

1  -1000.000  2.500  3.132  0.000  27.000  6.504  0.001 


38 


A  Numerical  Program  for  Flow  in  Eruptive  Conduits 


2 

-798 .867 

2.500 

3.132 

0.000 

21.601 

6.506 

0.001 

3 

-612.370 

2.500 

3.132 

0.000 

16.597 

6.509 

0.001 

4 

-459.127 

2.500 

3.131 

0.000 

12 . 486 

6.510 

0.001 

5 

-315.245 

2.500 

3.131 

0.000 

8.626 

6.512 

0.001 

257 

-0.082 

2.500 

8.959 

0.799 

0 . 447 

32.234 

1.000 

The  variables  are  listed  in  this  table  as  specified  in  the  input  file.  The  second  line  gives 
the  number  of  data  points  written  out  (257),  followed  by  257  lines  for  depths  extending 
from  the  base  to  the  top  of  the  conduit. 

Example  using  option  2  (specifying  pressure  gradient) 

For  the  second  example,  we’ve  taken  the  sample  input  file  and  changed  it  slightly  so 
that  the  pressure  gradient  (“pgrd”)  is  specified,  and  the  conduit  radius  (instead  of  log 
viscosity)  is  written  out  to  column  4  of  the  output  file.  After  typing  Conflow  to  start  the 
program,  the  program  echoes  the  input  variables,  then  prints  the  following  lines  for  the 
model  run: 


STARTING  RUN  NUMBER  1: 


mass  flux=  0.5244E+05  kg/s 


i  z  (m)  radius  log  Re 

1  -1000.000  2.500  2.319 

52  0.000  1.521  8.558 


vfgas  p  (MPa)  vel  (m/s)  mach  # 
0.000  27.000  1.000  0.000 

0.946  0.101  50.173  1.839 


successful  completion 


AFTER  ISENTROPIC  EQUILIBRATION  TO  1  ATM  PRESSURE: 
final  temperature  =  1143.42  deg.  C 
temperature  change  =  0.000  deg.  K 

enthalpy  change  =  0.0000E+00  J/kg 

max.  theoretical  velocity  =  50.17  m/s 


maximum  water  table  depth  that  will  allow  g.w.  influx  =  0.00  meters 

Negative  values  are  below  the  ground  surface, 
positive  values  are  above. 


Note  that  the  program  required  only  a  single  run  up  the  conduit,  because  the  p=pfinal 
boundary  condition  is  automatically  satisfied  at  the  final  depth.  Note  also  that  the  Mach 
number  (1.839)  at  the  conduit  exit  is  much  greater  than  1,  as  it  can  be  with  a  variable 
conduit  geometry.  Similarly,  the  exit  velocity  (50.17  m/s)  is  equal  to  the  maximum 
theoretical  velocity  because  the  erupting  mixture  has  fully  equilibrated  with  atmospheric 
pressure  by  the  time  it  reaches  the  surface.  The  maximum  water-table  depth  that  will 
allow  groundwater  influx  is  zero,  because  the  surface  is  the  only  place  where  the  two 
pressure  curves  (the  hydrostat  and  the  conduit  pressure  curve)  intersect. 


Closing  Comments 


This  report  is  intended  to  give  a  concise  summary  of  the  underlying  principles  of 
this  program  and  of  its  potential  applications.  It  will  probably  evolve  with  time  into 
something  more  complicated  and,  hopefully,  more  realistic.  If  you  intend  to  make 
extensive  use  of  this  program  or  would  like  to  find  out  about  new  revisions,  you  are 


encouraged  to  contact  the  author  at  (360)  993-8925  (e-mail  at  |1  u  mast  in  @  us  gs .  gov) 
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Appendix  A:  Calculating  M  eltThermodynamics 

Following  the  methodology  of  Ghiorso  and  Sack  (1995),  a  naturally  occurring  melt 
of  a  given  composition  can  be  considered  to  be  a  combination  of  certain  end-member 
components,  listed  in  Table  Al.  For  each  non-aqueous  component,  c  (the  specific  heat 

in  J/mole)  of  the  melt  has  been  determined  experimentally  and  is  assumed  to  be 
independent  of  temperature  (detailed  studies  of  individual  melts  (Neuville  et  al.,  1993) 
show  a  weak  temperature  dependence).  From  experimental  data,  Ghiorso  and  Sack  also 
give  the  molar  volume  of  each  component,  v  ,  and  its  derivatives,  (dv  /  dT )  ,  (dv  /dp) , 

(d2v  /dTdp) ,  and  (d2v/dp2)  at  conditions  of  T=  1673  K,  p=l  atm. 


Specific  heat,  density,  and  coefficient  of  thermal  expansion  of  the  melt 

Ghiorso  and  Sack  (1995)  use  a  regular  solution  model,  which  assumes  that  the  molar 
heat  ( c  ),  volume  ( vm ),  and  partial  of  volume  with  temperature  of  the  melt  are  weighted 

sums  of  the  properties  for  component  i: 


c  =  x  c hq 

pm  i  pi 

i= 1 
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i-w" 


i=l 


<K 

dT 

dp 


(  ^77  O  \ 


1=1 


i= 1 


dv 


v3G 

1  ^ ' 

Kdp  J 


(Al) 

(A2) 
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where ; q  is  the  mole  fraction  of  each  component.  The  overlines  indicate  that  these  terms 
are  given  per  mole  rather  than  per  kilogram  of  melt. 


Table  Al :  Names  and  formulas  of  end-member  components  in  melt 

Formula  of  end-  Component  name 

member  component _ 


Si02 

amorphous  silica 

Ti02 

rutile 

ai202 

corundum 

Fe203 

hematite 

Fe2Si04 

fayalite 

Mg2Si04 

forsterite 

CaSi03 

pseudowollastonite 

Na2Si03 

sodium  metasilicate 

KalSi04 

kalsilite 

H20 

water 
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Conflow  converts  these  terms  from  per-mole  values  to  per-kilogram  values  by 
multiplying  each  by  the  average  molecular  weight  of  the  melt  ( Mm ),  in  kg/mole: 

=  X>,M,  (A5) 

i= 1 
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pm 
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coefficient  of  thermal  expansion  = 


1  <K 

vmMtot  dT 


(A8) 


bulk  modulus,  Km- 


M 


tot 


Pm(dvm/dp) 


(A9) 


where  M ;  is  the  total  mass  of  component  i  per  mole  of  melt. 

Thermodynamic  properties  ofthe  melt 

Ghiorso  and  Sack  (1995),  use  a  regular  solution  model  to  calculate  the  main 
thermodynamic  properties  of  silicate  melts.  Their  formula  for  the  molar  Gibbs  free 
energy  of  the  melt  follows: 

n  n  n  n 

g  =  +rTYj  A  ln  A  +  - L E  w!7 A xj  + 

i= 1  i= 1  "  i= 1  7=1 

RT[xw  ln  xw  +  (1  -  .rH,)  ln(l  -  xw )]  (A10) 

The  subscripts  i  and  j  represent  individual  end-member  components.  The  variable  Wij 
represents  an  interaction  coefficient  between  components  i  and  j  in  the  melt.  Ghiorso  and 
Sack  (1995)  have  estimated  the  values  of  these  interaction  coefficients,  and  we  use  their 
values  in  my  calculations. 

Using  the  identity  s  =  -(dg  IdT )  (Moran  and  Shapiro,  1994,  p.  480),  the  previous 
equation  can  be  differentiated  to  obtain  the  molar  entropy: 

s=YJXis°-RYjxi\nxi  -R[xwlnxw  +  (1  - .rjln(l - *j]  (All) 

/=1  i= 1 
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And,  using  the  identity  h  =  g  +  Ts  ,  we  can  calculate  h  : 


i= 1  ^  i=l  j= 1 


(A  12) 


These  formulas  require  numerical  values  of  the  molar  enthalpy,  entropy,  and  Gibbs  free 
energy  of  each  end-member  component,  expressions  for  which  are  given  below. 


M  olar  enthalpy  of  each  component 

If  the  enthalpy  of  each  component  at  the  reference  temperature  (  7,-298  K)  and 
pressure  (p,= 1  atm)  are  known,  the  enthalpy  at  other  values  of  7  and  p  can  be  calculated 
from  the  following  formula: 


h°  =h°  + 

nT  p  nTr.pr  ~ 
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Substituting  (dh  /  c)T )  p  =  c p  and  (dh  / dp)  ,  =  v  -  (dh  /dT)p  (Moran  and  Shapiro,  1994,  p. 
491),  we  get: 


1  fusion 

K„.r.„ +  f^'dT  +  T^  + 
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(A  14) 


where  c  s°l  and  c1’'1  are  specific  heats  of  the  solid  and  liquid  phases,  respectively; 

A SfuSion,p  's  the  molar  entropy  of  melting  at  1  atm  pressure,  and  Tfusion  is  the  melting 

temperature  at  pr.  The  first  term  on  the  right  represents  the  enthalpy  of  formation  the 
component  from  the  elements  at  p=  1  atm,  7=298  K.  The  second  term  represents  the 
increase  in  enthalpy  as  the  mineral  is  heated  to  its  melting  point  at  1  atm  pressure.  The 
third  represents  the  enthalpy  of  fusion  of  the  component  at  1  atm  pressure,  and  the  fourth 
represents  the  enthalpy  added  to  the  melt  by  heating  from  the  melting  point  to  the  final 
temperature. 

The  terms  to  the  right  of  the  final  integral  in  Eq.  (A14)  represent  the  enthalpy 
change  associated  with  increasing  the  pressure  of  the  melt,  at  the  final  temperature,  from 
1  atm  to  the  given  pressure.  It  is  calculated  by  dividing  it  into  its  two  components:  (1)  the 
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integral  of  vdp ,  and  (2)  the  integral  of  -  (dv  /  dT)T  dp  .  The  term  vdp  is  integrated  as 
follows  (Ghiorso  and  Sack,  1995): 
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where  the  molar  volume  of  each  component  ( v 


1673  K,p, 


)  and  each  of  its  partial  derivatives 


have  been  tabulated  (Ghiorso  and  Sack,  1995).  The  derivatives  of  molar  volume  are 
determined  experimentally  and  are  approximated  as  independent  of  temperature  and 
pressure.  This  approximation  makes  it  possible  to  solve  the  second  term  as 


Pf 

-  J (dv  /  d p)Tdp  =  -(dv  /  d p)T  ( pf  -  pr )  (A16) 

Pr 


Because  the  melt  is  in  a  liquid  state,  each  component  is  assumed  to  be  a  liquid 
regardless  of  whether  the  melt  temperature  lies  is  above  or  below  the  fusion  temperature 
for  that  particular  component.  The  resulting  enthalpy  of  each  component  is  termed  an 
apparent  enthalpy  by  Ghiorso  and  Sack  (1995),  reflecting  the  fact  that  such  end-member 
components  may  not  necessarily  exist  as  separate  melts  at  the  temperature  and  pressure 
specified. 


M  olar  entropy  of  each  component 

Like  molar  enthalpy,  the  molar  entropy  can  be  calculated  by  integrating  from  the 
standard  state  of  T= 298  K,  p=\  atm: 
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where  s^  is  the  absolute  entropy  of  the  component  at  the  reference  temperature  and 


pressure.  Substituting  (ds°  /dT)p 
Shapiro,  1994,  p.  491),  we  get: 


cp/T  and  ( d.s  / dp)7 
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The  values  of  s°r  ,csp°l  and  cpq  are  obtained  by  experiment  (Berman,  1988;  Ghiorso  and 
Sack,  1995).  The  value  of  (dv 0  / dT ) p  is  also  obtained  by  experiment  (Ghiorso  and  Sack, 
1995),  and  is  approximated  as  being  independent  of  temperature. 

M  olar  Gibbs  free  energy  of  each  component 

The  molar  Gibbs  free  energy  (g°  )  is  calculated  from  the  molar  enthalpy  and 
temperature  using  the  definition  of  g  : 
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Values  of  g°  were  calculated  at  the  given  temperature  and  pressure  using  the 
following  formula,  from  Nicholls  (1980): 
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where  A  and  B  are  constants  whose  values  were  determined  from  experimental  data 
by  Nicholls.  The  value  P?  p  is  the  integral  of  specific  volume  with  pressure: 
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where  (Following  Burnham  and  Davis,  1974)  p  is  in  bars,  and  T  is  in  Kelvin.  The  above 
equation  uses  partial  molar  volumes  of  water  in  melt  that  were  determined  experimentally 
by  Burnham  and  Davis  (1971,  1974)  for  albite.  For  T=973-1473  K  (700° -1200°  C)  and 
p=l-1000  bars,  results  from  Burnham  and  Davis  give  v°wm  =20-35  cm  /mole — somewhat 
more  than  the  18  cm  /mole  for  liquid  water  at  25  °  C.  More  recent  studies  (Ochs  and 
Lange,  1999)  give  somewhat  smaller  values  of  vw  (-15-25  cm3/mole). 

The  MELTS  program  of  Ghiorso  and  Sack  (1995)  uses  values  of  A  (-33676.0 
J/(mole  K))  and  B  (18.3527  J/mole))  that  are  optimized  for  data  on  water  solubility.  The 
above  equation  indicates  that  the  chemical  potential  of  water  in  the  melt  is  a  linear 
function  of  absolute  temperature;  and  because  (dpi  IdT)  =  s*w ,  this  relationship  implies 
that  the  molar  entropy  of  dissolved  water  does  not  change  with  pressure  at  p= 0. 
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To  calculate  the  chemical  potential  of  water,  we  modify  the  equation  of  Nicholls, 
above,  by  subtracting  out  the  Gibbs  free  energy  of  supercritical  water  tabulated  by  Robie 
et  al.  (1978)  and  adding  the  identical  quantity  from  Haar  et  al.  (1984).  This  insures  the 
thermodynamic  properties  of  the  1TO  component  are  internally  consistent  with  those  of 
the  supercritical  fluid. 

The  enthalpy  and  entropy  of  dissolved  water  are  calculated  from  the  following 
thermodynamic  relations: 
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h  °  =  e°  +  Ts° 
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Partial  M  olar  Properties 

The  partial  molar  enthalpy  of  each  component  i  can  be  estimated  as  follows: 
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The  partial  molar  entropy  of  each  component  except  H2O  is: 


5,.  =s?  -Rlnxi  -Rln(l-xw) 


(A25) 


For  H2O,  it  is: 
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The  partial  molar  Gibbs  free  energy  of  each  component  except  H2O  is: 

n  ^  n  n 

Pi  =P°  +RT\ lnx,  +  ln(l  - .i'lv )]+  E wuxi -aEE^^./^ 

j= 1  Z  j= 1  k= 1 


For  water  dissolved  in  the  melt,  it  is 

Pw  =  PI  +  2R T In xw  +  E wwj x]  -^EE ^ (A28) 

2=1  z  7=1  t=i 

Equation  (A28)  is  used  to  calculate  gas  solubility  in  the  melt.  When  the  melt  is  saturated 
with  gas,  this  term  should  equal  the  Gibbs  Free  Energy  of  pure  H2O  gas  at  the  same 
pressure  and  temperature. 
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Appendix  B:  Calculating  the  Capillary  Number 

The  capillary  number  in  a  cylindrical  conduit  can  be  calculated  as  long  as  the 
strain  rate,  the  average  bubble  radius,  the  melt  viscosity,  and  the  melt  surface  tension  can 
be  evaluated.  The  melt  surface  tension  is  taken  to  be  0.34  N/m  (Proussevitch  and 
Sahagian,  1996).  Melt  viscosity  is  calculated  as  described  in  the  body  of  this  report. 

Other  values  are  calculated  as  follows: 

Shear  strain  rate.  For  laminar  flow  in  a  cylindrical  conduit,  the  velocity  profile  is 
parabolic  and  follows  the  relation: 
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where  u  is  the  velocity  at  a  given  distance  r  from  the  conduit  center,  and  R  is  the  conduit 
radius.  The  average  velocity  in  the  conduit,  u,  is  equal  to  wmax/2  (Bird  et  al.,  1960,  p.  46). 
The  shear-strain  rate  at  a  given  distance  r  from  the  conduit  center  is  equal  to  durldr : 
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The  average  shear  strain  rate,  £rz ,  is  equal  to  the  shear  rate  e'l ,  integrated  from  r= 0  to 
r=R  over  cylindrical  shells  of  infinitesimal  cross-sectional  area  27rrdr,  divided  by  the 
cross  sectional  area  of  the  conduit: 
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Integrating  this  equation,  the  absolute  value  of  the  average  shear  strain  is  £  =2u/3R 
Average  bubble  radius.  The  average  bubble  radius  (r)  is  calculated  from  the 
volume  of  gas  per  unit  volume  of  melt  ( vg/vm ),  divided  by  the  number  of  bubbles  per 
cubic  meter  of  melt  ( N \  the  "bubble-number  density").  The  quotient  (Vg/(Nvm)),  which  is 
the  average  volume  per  bubble,  is  converted  to  r  using  the  formula  for  the  volume  of  a 
spherical  bubble.  For  basaltic  lava-fountain  tephra,  Mangan  and  Cashman  (1996)  have 
measured  bubble-number  densities  of  about  1010  bubbles  per  cubic  meter  of  melt. 
Cashman  et  al.  (2000)  report  number  densities  of  1014-1016/m3  for  tephra  clasts  from 
silicic  plinian  eruptions.  Therefore  we  calculate  an  approximate  bubble-number  density 
as  follows: 
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where  mSiQ  is  the  mass  fraction  SiC>2  in  the  melt. 
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Appendix  C:  M  odifying  the  Source  Code  to  use  Papale's  Fragmentation 
Criterion 


To  use  Papale’s  fragmentation  model,  you  will  need  to  obtain  the  Fortran  source 
code  to  Conflow,  which  can  be  obtained  from  a  subdirectory  of  the  ftp  site 
“Elektra.wr.usgs.gov/Ftp_Access/Pub/Conflow.  Edit  the  source  code  as  follows: 

1)  In  the  main  program,  named  “main.f ’,  on  or  around  line  111,  remove  the  “c”  from 
the  beginning  of  the  following  line  of  code: 

c  read  (8,*)  ifragtype 

2)  In  the  same  file,  on  or  around  line  140,  add  a  “c”  to  the  first  column  of  the 
following  line  of  code: 

ifragtype  =  2  (sets  fragmentation  type  to  vfgas=0.75 

The  lines  of  code  that  determine  whether  the  melt  has  fragmented  are  on  or  around  lines 
703-710,  and  look  like: 

C  CHECK  TO  SEE  IF  WE'VE  REACHED  THE  FRAGMENTATION  DEPTH  YET. 

select  case  (ifragtype) 
case  (1) 

smod  =  2.5d+10  (elastic  modulus 

if  (dvdz . gt .  ( 0 . 01*smod/eta) )  ifrag=l 
case  (2) 

if  (vf gas . gt . 0 . 75 )  ifrag=l 

end  select 


The  variable  smod  refers  to  the  elastic  modulus,  eta  refers  to  the  bulk  viscosity,  and  dvdz 
is  the  extensional  strain  rate  in  the  conduit.  You  may  wish  to  judiciously  consider  their 
values  and  how  they  are  calculated. 

3)  To  the  input  file  conin,  Add  the  following  line  immediately  below  the  line  that 
gives  the  input  value  for  wt%  K20: 

1  fragmentation  criterion  (=1  for  Papale,  2  for  vfgas=0.75) 

To  use  Papale’s  fragmentation  criterion,  make  sure  that  the  integer  at  the  beginning  of 
this  line  equals  1.  To  use  the  ty,=0.75  fragmentation  criterion,  change  it  to  2. 

4)  Recompile  the  program.  A  readme .  txt  file  in  the  ftp  directory  containing  the 
source  code  will  tell  you  how  to  do  this.  You  will  now  have  the  option  of  choosing  either 
Papale’s  fragmentation  criterion  (by  entering  a  “1”  at  that  line  of  the  input  file),  or 
vg=0.75.  However,  you  will  not  be  able  to  use  this  version  with  the  Windows®-based 
front  end.  You  will  have  to  run  your  models  from  the  DOS  window  (if  using  a  DOS  or 
Windows  -based  machine)  or  from  the  command  line  (for  other  operating  systems. 

Appendix  D:  Calculating  umax 

If  the  erupting  mixture  reaches  the  conduit  exit  before  the  pressure  has  dropped  to  1 
atmosphere,  it  will  abruptly  expand  to  equilibrate  with  atmospheric  pressure.  Expansion 
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and  adiabatic  cooling  will  accompany  this  decompression.  If  these  processes  take  place 
without  frictional  dissipation  of  energy,  the  process  is  said  to  be  isentropic,  and  maximum 
amounts  of  acceleration  and  cooling  can  be  calculated.  In  this  program  the  calculations 
are  done  with  the  assumption  that  the  mixture  acts  as  an  ideal  “pseudogas”  (Kieffer, 

1984).  That  is,  the  mixture’s  bulk  properties  approximately  follow  the  ideal  gas 
relationship,  pv=nRT.  For  ideal  gases  and  pseudogases  expanding  under  adiabatic, 
isentropic  conditions,  the  pressure  and  temperature  before  and  after  decompression  are 
related  by  the  equation  (Moran  and  Shapiro,  1992,  p.  104): 
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where  the  subscript  f  refers  to  the  “final”  value  in  the  conduit,  before  decompression,  and 
“e”  refers  to  the  value  after  decompression.  Temperatures  are  absolute.  The  variable  y  is 
the  ratio  c/cv,  where  cp  and  cv  are  the  specific  heats  at  constant  pressure  and  constant 
volume,  respectively,  of  the  magma-gas  mixture.  Those  specific  heats  are  given  by  the 
equations: 
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where  cPig,  cVig ,  are  the  specific  heats  at  constant  pressure  and  volume,  respectively,  of  the 
gas  phase,  cpjn  and  cv,m  are  the  specific  heats  at  constant  p  and  T,  respectively,  of  the 
liquid  magma,  and  cPiX  and  cVyX  are  the  analogous  values  for  the  crystal  phase,  we  assume 
that  cPiin=cv,m  and  cPtX=cv,x.  Specific  heats  of  the  gas  phase  are  calculated  using  the  method 
of  Haar  et  al.  (1984).  The  specific  heat  of  the  liquid  magma  is  calculated  using  the 
methodology  of  Ghiorso  and  Sack  (1995;  Appendix  A). 

Once  the  adiabatic  temperature  change  has  been  calculated,  the  change  in  specific 
enthalpy  (h)  of  the  mixture  during  decompression  is  computed  from  the  following 
equation  for  ideal  gases  (Moran  and  Shapiro,  1992,  p.  96): 


he-hf=  cp(Te  -7f)  (D4) 

In  addition  to  assuming  ideal  gas  behavior,  this  equation  assumes  that  cp  is  invariant  over 
the  range  of  temperatures  experienced  during  decompression. 

The  maximum  theoretical  velocity  is  then  calculated  assuming  that  all  of  the  change 
in  enthalpy  of  the  expanding  mixture  is  transformed  into  kinetic  energy.  This  implies  that 
an  insignificant  amount  of  energy  goes  into  lifting  of  the  material  or  to  frictional  heating. 
Given  these  assumptions,  the  maximum  theoretical  velocity  (wmax)  is: 

Umax  =  V U)  +2(hf  ~he)  (D5) 
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